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ABSTRACT 


The axtsyminetric heat equation. resulting from a point-souce of heat applied 
to a metal bloek. 1s solved numerically; both iterative and multilevel solutions are 
computed in order to compare the two processes. The continuum problem is dis- 
cretized in two stages: fimte differences are used to discretize the time derivatives. 
resulting moa fully impheit backward tnne-stepping scheme. and the Finite Volume 
element (IVE) method is used to discretize the spatial derivatives. The application 
of the EVE method to a problem in cylindrical coordinates is new. and results mn 
stencils which are analyzed extensively. Several iteration schemes are considered, in- 
Cluding both Jaeobi and Gauss-Seidel: a thorough analysis of these schemes ts done, 
using both the spectral radi of the iteration matrices and local mode analysis. Us- 
ing this discretization, a Gauss-Seidel relaxation scheme is used to solve the heat 
equation iteratively. A multilevel solution process is then constructed. including the 
development of intergrid transfer and coarse grid operators. Local mode analysis 1s 
performed on the components of the amplification matrix, resulting in the two-level 
convergence factors for various combinations of the operators. The multilevel solu- 
Lion process 1s implemented by using multigrid V-cycles; the iterative and multilevel 
results are compared and discussed in detatl. The computational savings resulting 


from the multilevel process are then discussed. 





DISCLAIMER 


Phe computer program in Appendix B is supplied on an “as ts” basis. with no 
warrantees of any kind. The author bears no responsibility for anv consequences of 


using this program. 
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i INTRODUCTION 


‘\ topic of general interest is the numerical solution of partial differential equa- 
tious (PDEs). such as the Navier-Stokes equation (Equation II.1). Many aspects of 
Iliese types of problems present interesting challenges to constructing numerical so- 
litions., sueh as non-linearities, high-speed flows that cause numerical peculiarities, 
unusnal boundary conditions. uncommon geometries, and so on. Developing pro- 
cesses that solve these tvpes of problems is difficult. and often results in complicated 
solutions schemes which are costly to implement. .\ common goal is to attempt to 
streatutine existing solution processes, or develop new ones. in order to reduce com- 
putational complexity and cost. 

At issue is how to transform a continuum problem into a discrete one (dis- 
cretization), and how to construct a numerical process that will solve the discrete 
problem. The finite volume element (FVE) method (see [Ref. 1]) has proven to be 
a useful tool in developing discretizations, particularly for problems that require the 
enforcement of conservation laws. One of the more successful methods for streamlin- 
ing the solution of PDEs, particularly elliptic equations, is a multilevel technique, to 
wit. multigrid (see (Ref. 2], [Ref. 3], [Ref. 1]). While this method enjoys remarkable 
success in solving certain classes of both linear and non-linear PDEs, its applicability 
to solving other types of equations awaits development. The goal of this work 1s two- 
fold: to apply the FVE method to discretize a particular equation, and to implement 
multigrid in solving the resulting discretized problem. 

The approach that is taken is to begin with a specific physical problem which 
results in the Navier-Stokes equation, and consider a subsidiary problem, namely the 
heat equation; due to the nature of the problem, cylindrical coordinates are used. 
The resulting problem, to which the FVE discretization and multilevel solution are 


applied, is the axisymmetric heat equation, involving the use of a point source. More 


specifically, finite differences are used to discretize the temporal portion of the prob- 
lem, resulting in the use of a fully implicit backward time-stepping scheme, and the 
Wl method is used to discretize the spatial part. The application of the FVE 
inethod to a problem in cylindrical coordinates is new. and results in stencils which 
are analyzed extensively. Once the discretization has been constructed, a number of 
relaxation schemes. including both Jacobi and Gauss-Seidel, are considered; a thor- 
ough analysis of these schemes is done, using both the spectral radii of the iteration 
matrices and local mode analysis. The goal is to choose one of the relaxation schemes 
fur use in an iterative solution of the problem; Gauss-Seidel is the method of choice. 
The specifies of the multilevel technique are then developed and analyzed. Local 
inode analysis is performed on the components of the amplification matrix, resulting 
in the two-level convergence factors for various combinations of the operators. The 
multilevel sohition process is implemented by using multigrid V-cycles; the iterative 
and multilevel results are compared and discussed in detail. The computational sav- 
ings resulting from the multilevel process are then discussed. Insofar as multigrid 1s 
quite successful in a number of different contexts. the results of our work are some- 
what disappointing in that we are not yet able to achieve the same level of success. 
On a more positive note, thts work applies the FVE method to a problem in cylindrt- 
cal coordinates, and makes use of Maple software to compute the necessary integrals 
to construct the discretization (see Chapter III and Appendix B). Additionally, the 


techniques that are applied do result in a solution to the problem. 


le PROBLEM STATEMENT 


A. BACKGROUND 


Our ultimate goal is to be able to numerically solve the Navier-Stokes equation. 
Doe 2 
—+u-Vu= —--Vpt+uwvt. (IT.1) 
p 


where wis the velocity vector, ¢ is the time. p is the density, p is the pressure. and v 
is the kinematic viscosity. One approach to such solutions ts to consider this problem 
as a collection of smaller problems. By solving each of the sinaller problems, we 
can assemble the collection of solution processes to solve the original problem. The 
purpose of this work is to focus on an initial step in the eventual construction of a 
numerical solution to Equation II.1. 

One way to subdivide a problem such as this into smaller problems ts to 
consider a specific physical problem that gives rise to the Navier-Stokes equation, and 
then isolate the different aspects of that problem. To that end, we consider a semi- 
iutinite block of metal, with a heat source applied to the horizontal surface; above the 
horizontal surface ts an (invisctd) nonconducting gas. The result is a pool of molten 
metal with a free surface, surrounded by the solid portion of the unmelted block (see 
ieure 1). The total heat flux Q is constant, and far away the solid approaches the 
uniform cold temperature, T,. The resulting thermal and flow fields are assumed to 
be axisymmetric and steady, and are governed by conservation of energy in the sold 
and by conservation of energy, momentum, and mass in the pool. We therefore have 


the following system of equations: 


IT | 

solid: a. = «VT (II.2) 
i | 

liquid —: a +i-VT =nV’°T ORS 

= +i. Va = —-Vp40V7a (11.4) 

er (11.5) 


where Tis the temperature. « is the thermal diffusivity, and wu. t, p, p, and v are as 


above, 





solid 


Figure |. The molten-pool problem. 


The portion of the above problem which is our focus is the conduction of heat 
in the solid, governed by the heat equation (Equation IJ.2), to which the following 


boundary conditions apply: 


sumace a> = (hee Ke = —q(r) (11.6) 
axis.) f) = 0uee gh = (11.7) 
Or 
far away ryz—7o0oo : Tok, (11.8) 


where & is the thermal conductivity, and g(r) is the imposed surface heat flux (large 
at r = 0, falling off to zero at some small value of r, such that {5° q(r)2ardr = Q); the 
condition in Equation II.7 results from the axisymmetric nature of the problem. Since 
the ultimate goal is to solve the Navier-Stokes equation in cylindrical coordinates, the 


heat equation will also be solved in this coordinate system.(see [Ref. 4]) 


B. POINT-SOURCE PROBLEM 


Our goal is to take the first step toward solving the Navier-Stokes equation 
(Iequation IT.1]) that arises from considering the molten-pool problem by solving the 
associated heat equation (Equation 1.2). We therefore assume cylindrical geome- 
try and azimuthal svmmetry. We further simplify the conditions imposed on Equa- 
tion IL.2 by considering that the heat source applied to the horizontal surface of the 
block is a point source, with the total heat flux Q = 1. and that far away the solid 
approaches the wmform cold temperature of 7. = 0. That is, the imposed heat flux 
me} = «(7 ) as the Dirac celta function. with J, qtr )j2ardr = 1. 

The existence of the solution to this type of problem is well established. How- 
ever, as Is the case in general when Neumann boundary conditions apply. absent a 
specified initial temperature distribution. the solution is not unique (see [Ref. 5]). 
While we are interested in solving the point-source problem numerically in cvlindrical 
coordinates, it has spherical symmetry, so it can be solved analytically in spherical 
coordinates more easily. Therefore, for the analytic solution, we rewrite Equation I[.2 
in spherically symmetric coordinates: 

OT ieee, Oia 
ees (IT.9) 
a an (aR) 
where ? is in spherical coordinates, R = Vr? + 2%. Beginning with an initial tem- 
perature distribution of 7’ = 0 everywhere and Q = 1, the solution to (11.9) with 


boundary conditions 








Siibliace pace — aU) © a = —q(r) = —d(r) CUE) 
z 
axis) te Ue = ae (11.11) 
Or 
far away r,z7oco : T7>0, (ii) 
is found to be 
td) = f a Uses) 
st) = 5 perte Weak 


where erfe(x) is the complementary error function. 
2 eee 
eric(x =) —ent(<),. fen — | eds. 
0 


A\s t + oo, the steady solution becomes 


l 


Orin 


T(R) (1.14) 


(for a derivation of this result, see [Ref. 5]). In (1, z) coordinates, where R = Vr? + 2?, 


these become 


Z 2 

TG a Le (11.15) 
InKVre + 2 D7 ee 

i) eae (11.16) 


for computational purposes, the idealized problem of a semi-infinite solid 
is truncated to a finite domain in cylindrical coordinates with azimuthal symme- 
try. At the far boundaries on this finite domain, homogeneous Dirichlet conditions 
are imposed to approximate the conditions in the unbounded problem. To further 
simplify computational work, we assume that the (r,z) domain is the unit grid 
Q = [{0,1]x{0,1] € R*;we will eventually assume x = 1. Thus, the equation we 


wish to solve is 


OT | 
— =KV°T, Uy 
Re ahi 

with boundary conditions 
surface z=0 aa = —d(r) (11.18) 
axis, 7 — gh a) (IT.19) 
Or 

far boundaries 332 =) a0 eee (II.20) 


ie DISCRETIZATION 


ln order to solve the conduction problem numerically (Equation I[.2. with 
mounldary conditions Equations [1.18. 11.19. and 1.20). the equation representing the 
continuum problem must be discretized. That is. the values of the unknown function 
f are to be determined from a set of equations which in some sense approximate 
Squation [[.2. Ultimately, a discrete approximation to the heat equation should meet 


the following criterta: 


1) Provide for a unique solution to the problem. 


2) The solution should be “close” to the exact solution for “sufficiently small” 
vyrid spacings. 


3) The solution should be “effectively computable.” 


The significance of property |) is obvious; property 2) relates to the questions of 
convergence and consistency for the discretization scheme. Property 3) relates to: 
a) the amount of computational work required to solve the problem, and b) the 
behavior of roundoff errors in the computed solution. The growth of the roundoff 
error Is related to the notion of the stability of the discretization scheme. (see [Ref. 
i) 

Solution of the heat equation (Equation II.2) requires treatment of both space 
and time derivatives. The Finite Volume Element (FVE) method is used to discretize 
the spatial portion of the problem; finite differences are used to discretize the tem- 
poral portion. The approximation properties of the FVE method are fundamentally 
different from those associated with the finite difference method. Hence, this ap- 
proach treats time and space differently. The goal in applying the FVE method to 
produce spatial discretization is to make use of the advantages of the Finite Volume 
(FV) method, its ability to be faithful to the physics in general and conservation 
in particular, coupled with the flexibility of the finite element representation of the 


unknown functions (see [Ref. 1]). As outlined in [Ref. 1], the basic approach is 


io partition the spatial grid in two wavs: as the union of a set of finite elements. 
the vertices of which comprise the grid points on which the unknowns are defined 
(see Figure 2), and as the union of a set of control volumes (see Figure 3), one for 
each grid point (the boundary points require separate treatment. as will be discussed 
later). The elements are used to form basis functions, a linear combination of which ts 
used to approximate the unknown function. Upon substitution of the approximation 
ito the equation to be solved. integrals over each control volume are taken. The 
integrals enforce conservation on each control volume, and therefore the partition 


enforces conservation on the entire domaim. The system of equations generated by 


iutegration yields the discretization of the problem. 


VAVAVAVA 


NIN 
NIA 


BVAVAVAVA 


2 N 


P __—______._» 


Figure 2. Partitioning of problem domain cross-section 22 by elements. 


The process of discretizing Equation I1.2 in both space and time is complicated, | 
the application of the FVE method to the space derivatives being the more complex 


portion of this process. Therefore, before proceeding to a discussion of the PVE 
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Figure 3. Partitioning of problem domain cross-section 2 into sub-volumes. 


method, the finite difference method is applied to the time derivative. In order to 
facilitate this presentation, we present a one-dimensional example of this type of 


discretization (see {Ref. 7]). Consider the one-dimensional version of Equation II.2. 


a GPUE 


eee IIf.1 
Us a eof ves 


where « = 1. One finite-difference approximation to Equation III.1 ts 


ee en ae ee ee ee 
g h? 


where T is the exact solution of the approximating difference equations, 
2 el 1 mrs Uae cape eat) Cem or (saa ne ee 
Let r = 4, and this can be written as 


jae = a ae = (1 a aes = FL ey as 


Thus. we can compute the unknown T7;,.,4; at the (nr + 1)st time step using the 
known values from the nth time step, giving an explicit formula for determing the 
unknowns.|Ref. 7] 

While this explicit relation provides a simple means to compute unknown 
values. 1s has a sertous drawback. The process is only valid for 0 < a < s, or 
es i and therefore the time step A? = g 1s necessarily small. Moreover, in order 
to achieve reasonable accuracy, Ax = h must also be kept small. There is, however. 
a method which reduces the total amount of calculation and is valid (1.e., convergent 
aud stable) for all finite values of r, which was proposed by Crank and Nicholson 
(see (Ref. 7]). The technique is to consider Equation III.1 as being satisfied at 
the midpoint {Ah.(n + 4)g} and replace ~t by the average of its finite-difference 


approximations at the nth and (n+ 1)st time steps. That is. the equation 


ar) _ (aT 
Ot pee — Ox? a 


Lee ie. | ‘ay eee = aA re = p Airey ore cae os Zoe == lee 
Se 


is approximated by 
g 2 


which gives 
—7rT ei nti ar Z te Pa GS ace a a Vac) =, i Os i (2 ae 2) ae ae rT ep tyns 


where r = 74. Therefore, we now have three unknowns at the (n + 1)st time step 
written in terms of three known values at the mth time step. Thus, the Crank- 
Nicholson method generates a set of equations which must be solved simultaneously; 
it is an implicit method. 

Although the Crank-Nicholson method for Equation III.1 is stable for all pos- 
itive values of r in the sense that all errors eventually tend to zero as n tends to 
infinity, large values of r, e.g., 40, can introduce oscillations in the solution (see [Ref. 


7|). Therefore, a more general finite-difference approximation to Equation III.1 is 


10 


found bv using a weighted average of the finite-difference approximations of Lt at 
tle wth and (2 + L)st time steps, which is given by 
Peel — fk Oe Nae Ooo a ee eae 0 ee Oe er 
BEE )+(1-e) (Se) 
gy He fe 
where 0 <a <1 is possible. Note that a = 0 gives the explicit scheme. a = 4, Crank- 
Nicholson, and a@ = | a fully implicit backward time-difference method. The scheme 
is unconditionally valtd (stable and convergent) for 5 a7 <— |, but ior Ua = , we 
must have r <1/(2 —da). Thus. for example. for a = 0. r = oS 7 OF ig A is a 
requirement to guarantee validity.[Ref. 7] 
We now apply this fintte-difference method of discretizing the time derivative 
to Equation I1.2. Following the work done above. a general finite-difference approxi- 


mation for the time derivative of 


OT | 
— = KV’*T Pie? 
at ue) 
Is given by aes | | 
ie Tis ae + (I _ aed (LTS) 
g 


where g 1s the time step, and 0 < a < 1 allows for a weighted average of the current 
and future time steps. The current time step is designated by n, n + 1 ts the next 
time step; the value of a determines the type of time stepping. However, we have 
determined that in order to guarantee the validity (stability and convergence) of the 
expheit verston of this scheme (a = 0), we must have g < A We consider this 
value of g to be too sinall to be of computational use (which is verified by numerical 
experiments), and therefore we will use the values a = 4 and a = 1. Thus, we have 


the semi-discrete relationship between the current and subsequent time steps: 


(- = Ce ase = fe ais el = Gale ae (III.4) 
g 


This relationship is used to establish a discrete set of equations, which gives rise to 


the discretization of the problem. 


| 


A. BASIS FUNCTIONS: FINITE ELEMENTS 

We begin the discussion of discretizing Equation IT.2 in the spatial variables by 
recognizing that it is necessary to make an assumption about what types of functions 
mav be used to approxtmate the unknown function T; we consider that T can be 
approximated by continuous. piecewise linear functions. A basis ts then chosen for 
the space in which these functions are found. That is, a set of functions is selected 
whose linear combinations determine the space of functions of interest. Therefore. we 
choose basis functions. @,j;(7. 2), which are assumed to be continuous and piecewise 


linear. with 


eh! 
Hees) = > \ eee: (bles) 
k 


=O 7=) 

[he next step is to determine how to construct these basis functions; the element 
partition of the domain 22 is used as a framework for this purpose. (see Figure 2). 

Since we are assuming cylindrical geometry with azimuthal symmetry, the 
ltrecti 5 Ee ielees if id of izh=~(N= b 
directions of interest are r and z; a uniform grid of step size h = 3 (N = number 
of grid points) is applied in both directions. In this way, we can designate function 
values on the grid 2 by 


Th, = ie Sey, 


where 2, = kh, r; = jh, and k,j are the row and column indices. Each square 
is subdivided into two triangular elements by a diagonal (oriented in the direction 
of increasing 7 + z)'; the basis functions to be used in approximating the unknown 
function are constructed using these triangular elements. For each interior grid point, 
or node’, there are six associated triangular elements: NNE, ENE, SE, SSW,WSW, 
NW (see Figure 4). 


''The orientation of the elements is based on the necessity to apply this technique to solving the 
molten-pool problem. In particular, we will be required to track the movement of the phase-change 
boundary between solid and liquid; element orientation has been chosen to facilitate keeping track 
of this inoving boundary. 


*For a two-dimensional grid with N intervals in each dimension, there are (N - 12 interior grid 
points, (N + 1)? total nodes. 


1D, 


ik+by k+l y+) 


(kj-1) (k+l) 





Kel D (he Ty) 


BQ tT > 


eure 4. Hat function, top view. & is the row index (z direction), 7 is the column 
index (7 direction), y~ direction is into the plane of the paper (@). 


[f we suppose that this collection of elements ts joined along shared borders, 
and that the elements can be “stretched”, then the finite element for a particular 
node. say (k,7), is formed by anchoring the collection along its external boundary (in 
the (7, z) plane) and extruding the grid point in a direction perpendicular to the (r, z) 
plane and parallel to the tangent to the ~ direction at the grid point (A, 7). forming 


the famihar “hat” function (see Figure 5). Note that 


1 at the node (1r;, <4) 
Pk,j = 


0 at all other nodes 
and is piecewise linear over each triangular element, giving the hat shape. 
By choosing these hat functions to be our basis functions, the coefficients of 
the d,,(7, 2) in Equation III.5 become the values @,; = Ty,;. That is, we can now 
specify the nodal values of T’ on the grid as T,,; = T(r;, 2%), with 


N WN 
TGs 1 SS dE yeroyied Ce nal (IIT.6) 


k=0'97=0 


13 





Figure 5. Hat function, oblique view. 


Since the @,,, are piecewise linear, a continuous representation for function values of 


Tis found by linear interpolation of nodal values between nodes. For example, 


Up to now, we have considered the sampled values of the unknown 7’ at the 
(.V + 1)* grid points as an (N + 1)x(N + 1) two-dimensional array; it will later 
be convenient to consider 7 as a one dimensional array. One simple method to 


accomplish this is to order the elements T;,,; lexicographically. For example, 


To,0 To1 ss To,N+1 
Ti 0 Ti mi Ti N+1 
Tne+io Tein +++) Tn+i.n4i 


can be written as 
T 
(Too; Tan To,N41; Tor ieee LN aye Peal Ween Tn 41,1; ae Tn41,N+1) . 
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which we will normally consider to be a column vector. Thus, we can relabel the 


values ye, = fy, = (N+ lA + 74+ 1 and &.7 =0: V 4+ 1. 1n Equation IIT.6 so that 


(a Cees a — 3 iy te |, FO re 


Be FVE STENCILS 


With the construction of the basts fictions over the finite elements. we can 
how icorporate the control volumes to complete the finite volume element discretiza- 
tion. Below is an example the FVE method applied to a simpler problem: the purpose 
Is Lo motivate and explain the work that is necessary to handle the specifics of apply- 


ig the FVE method to Equation II.2. This explanation closely follows the example 


in [Ref. 1]. 


1. Example: 2—D Potential Flow 

The example problem we consider is the potential flow problem in Euclidean 
two-dimensional coordinates. The domain for the problem is the unit grid 2 = 
(0. 1]x(O.1] € R*, where x = 0 and y = 0 are the near boundaries and x = | and 


y = | are the far boundaries. The governing equation ts 
VatpV wb) =o, (III.S) 


where p is the density, y is the flow potential (i.e.. Y, and wy, are the components 
of the fluid velocity in the x and y directions), and 7 is the interior source flow rate. 


The boundary conditions are 


near boundaries : (pVw)-% = 7%, (Neumann boundary condition), 


far boundaries : w= wo (Dirichlet boundary condition), 


where n is the outward normal. Suppose that V is one of the control volumes in 2. 
similar to those depicted in Figure 3, where r and z are replaced by zx and y. Since 


this problem is in two dimensions, the control volumes V are actually areas, and the 
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corresponding surfaces 5’ are the perimeters of these areas. Integrating Equation IIL.8 


Over We aGe herve 


| V -(pVw)dV a | ndV. 
4 V 


By applying the Gauss Divergence Theorem, the volume integral on the left-hand 


side becomes a surface integral. 
| OV nds = | aa (111.9) 
S 


ln terms of the physical units associated with the potential flow, the relationship 1s 


mass length mass 
——— . ——— - area = ———————_ - volume; 
volume time time - volume 


each side represents a flow rate of mass per unit time, and (pVw)-n represents a flux 
across 5. Therefore, the net flow from the interior source in V is balanced by the 
flow across the surface S, which can be interpreted as a mass conservation law for the 
control volume V.[Ref. 1] 

By partitioning the domain 2 into a finite collection of control volumes (see 
Figure 3, where r and z are replaced by x and y), we can impose on each control 
volume the condition represented in Equation III.9 (the boundaries require special 
treatment: a discussion of these considerations is included later). As noted earlier, 
the imposition of conservation on each control volume results in conservation over 
the entire domain. This integral condition on each control volume will be used to 
compute the discretization; the number of discrete equations is the number of control 
volumes, say m, used to partition 2. To complete the discretization process, the 
unknown function w will be approximated by @ in R™. By appropriately choosing 
basis functions, we can construct an approximation v expressed in terms of its nodal 
values. Consider a triangular element partition similar to that in Figure 2 (where r 
and 2 are replaced by x and y) and let TY be the space of continuous, piecewise linear 
functions on 2 associated with this partition. Ignoring for the moment conditions at 


the boundaries, the FVE discretization of Equation III.8 comes from finding a v € T 
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seh that Equation TL.9 holds for each of the control volumes in the partition of Q. 


Phe problem im R™ ts then defined when e ts expressed in terms of nodal basis for T: 


rt 


i eS Toye (111.10) 


w= 

where uw, 1s the value of v at the wth node and ©,, is the “hat function” associated with 
the wth node (as above, we lexicographically order the nodes, resulting in a single 
subscript), Substituting Equation IIT.10 into Equation III.9, we have the matrix 
Cqnation 


byt (ITD) 


where £ 1s the nxn matrix with entries 
ia [ fevon ands (111.12) 
and 
fu = |) nav (111.13) 
Vw 
(except for boundary conditions).{Ref. 1] 

The treatment of the boundaries depends the type of boundary condition. 
Neumann conditions can be imposed indirectly by substituting the flux value 7; into 
the appropriate term in Equation III.9. Specifically, suppose we choose the control 
volume V that includes the origin (lower left-hand corner, Figure 3), where the surface 


vf the control volume coincides with the boundary along the two axes. The integral 


concition for V is 


/ ind + f SGD tee [ nav: (111.14) 
rv=0 ry HO 


Equation III.14 is substituted for the interior condition in Equation [11.9 into the 

discrete approximation v for the corner volume containing the origin.{Ref. 1] 
Dirichlet conditions are imposed directly; u, takes on the value of y, at each 

Dirichlet node, i.e., each node on the far boundaries (z,y = 1). This approach may 


lead to fewer unknowns u,, than equations, a problem easily resolved by discarding the 


ler 


cquations associated with the Dirichlet nodes. In this case. the collection of control 
volumes no longer partitions the domain and. therefore. conservation no longer 1s 
strictly enforced. Tlowever. it is at the Dirichlet nodes where the loss of conservation 
is nota concern in general.[Ref. 1] 

Assuming that there is a balance in the number of equations and unknowns, we 
now must implement a numerical rule for evaluating the integrals in Equations III.12 
and HE.13. For Equation [II.12. we use the following rule on each segment A that ts 


part of the interior surface S associated with a volume V: 


/ (COs ane / Ween 
A A 


where P4 is point of intersection of A and the grid Imes passing through the node 


associated with V. For Neumann boundary segments. we use the quadrature rule 


[vids = e(PadAl, 


where [A] is the “surface area”. i.e., length, of A. For Equation III.13, we use 


A ndV ~n( Ny)|V\, 


where Vy is the node associated with V and |V| = fy dV is the “volume”. 1.e., area, 
of V.{Ref. 1] 
Numerically evaluating Equations III.12 and III.13 yields the FVE discretiza- 


tion of Equation III.8. The value associated with each node, say (p,q), is a sum 


w=p+l1,z=q+l 


where the coefficients a,,, are determined by the integration. A convenient way to 


represent the sum associated with each grid point is to place the coefficients a,,, into 


a stencil, 


Ap+l,q 
Qpg-1 Ag ott | Upq = %p—1,gUp—1,g TF %pq—-1Up,q—1 


Op—1,g 


+ Ap gUyp,g + UptigUp+iyg + Up,g+1Upgtl 
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where, for this example. the corner values im the stencil are all zero. The value tor 


node tp.qg) is determmed by applying the stencil to u 


Pa 
To see what type of stencils are produced. consider the discretization on a 
auiforin mesh with grid size h = — in both coordinates. We use double subscripts 


pg =: m— 1, where 0 corresponds to the Neumann boundary nodes and m — | 
corresponds to the Dirichlet boundary nodes. Written m stencil form, the interior 
nodes. O < pq < m— |. for the equations tn Equation III.11 are as follows (3° 


indicates the sum of the outer stencil entries): 


p((p + 3)h, qh) 
Aphilq- +)h) -> p(ph.(q + 5h) Ung = h?n(ph. gh). 
p((p — 3). qh) 


lor the corner Neumann node (p,q = 0), the stencil is given by 


2 


h h h h 
—y>  4p(0, 2) | voo = = 140, 9) —5 (vs00 ha + 6s(G.0)) 


For both stencils, the coefficient of the unknown to which the stencil is applied is 
= a 

For grid points adjacent to a Dirichlet node, the stencil entry reaching to the 
Dirichlet node value is moved to the right-hand side as a coefficient of the boundary 
data. Otherwise, the equations for these points are the same as above. For example, 


at the pomt (0 < p< m—1,q=m-—1) we have 


p(ph, (q — 3)h) -> Upm—1 = h’n(ph,1 —h) 
p((p = sh, See h) 


h 
—p(ph, PS 5 Polph, Es 


where >> 1s the sum of the outer stencil entries without the boundary terms removed, 


a 


h 3, 1 1 
2d = p(ph, 1 — 5) + o(ph, (a - 5)h) + o((p + 5)hy1 —h) + a((p — 5), 1 — 4h). 
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(see Ref. 1]) 


Ne Application to Axisymmetric Heat Equation 

This technique is now used to compute the FVE stencils for Equation I[1.2. by 
combining the finite elements (see Figure 2) with the control volumes (see Figure 3). 
A control volume for the /th grid point. call it \y, is a toroidal prism (see Figure 6). 
[1 is eenerated by taking the two-dimensional sub-volume for that point (in the (7, z) 


plane). and sweeping through 360° in the direction. 


control 


volume 





Figure 6. Toroidal volume for the conduction problem. 


‘e 
Integrating Equation III.4 over all control volumes V;, where V = ae V 


partitions the domain 2, we have 


| — ~anV?|\T dV = / 1 4 (1 —a)eV?) T'aV, (111.15) 
Vv \g VAG 
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or, upon application of the Gauss Divergence Theorem the volume integral of the 
Y-T term becomes a surface integral. and we have 
~ | fy wi Pameae [ Sit ey ee - [ THlV + (1 —e) ef See SLA EOS 
u g 
where nis the outward normal. Substituting the expression for approximating the 


unknown function 7 (Equation [I].7). we have 


(N+1) l 
De * f grtdV — an f Vert. nds) (eiae ({1T.17) 
i=] ui ‘ S 
(N+1)° | 
= a “/ gdV + (1 —ajn | Mo i 
‘| oo ) 


where we now have integrals over each of the (V + 1)° control volumes, resulting in 
a set of (V+ 1)* equations. This set of equations can be written both as an operator 
equation, L[T"t'|] = f(T"), and as a matrix equation®?, M[T"t'] = f(T”), where the 


operator L and the matrix M are given by 


L(j-«") 


; / SAVIN, Mey [ Wounds) 


L 


M 


f(T") and Ges are given by 


es [.(G+0-e. Vv?) rev, 


i(P") » | [ stav + (1a) Jx | Var -nds| T. 


[=1 
Any function whose coefficients satisfy the resulting set of discrete equations neces- 
sarily satisfies the conservation law over any volume made up of the union of control 
volumes, except possibly at the boundaries. The Neumann conditions at r = 0 and 


z = 0 are incorporated indirectly into T, by substituting the (zero) normal derivative 


3The operator may be represented as a continuous operator, L; as a discrete operator, L’; or as 
a matrix, M. 


Zl 


value tnto the appropriate term in Equation [I[.15. The Dirichlet conditions at the 
lar boundaries are imposed directly on T: control volumes that would usually be as- 
soctated with these boundary nodes are elimiated (see Figure 7). Thus. the reduced 
collection of control volumes no longer partitions the grid Q. which slightly impairs 
conservation in the discretization. However, since we require the temperature to fall 


olf to zero at these points. this loss of conservation is not a concern.[Ref. 1] 


N 


to 





Figure 7. The problem domain cross-section is depicted after partitioning into sub- 
volumes. Homogeneous Dirichlet boundary conditions obviate the necessity to define 
volumes for grid points at the far boundaries. 


The integrals over the basis functions have been computed using Maple soft- 
ware and a program designed by David Canright (see Appendix B and [Ref. 8]). In 
order to calculate the integrals, we must be able to represent the unknown function 
over the basis functions, which themselves are collections of triangular planes. There- 


fore, the method is to determine first the function for the plane through three given 


[ie 


potts (the cases of vertical planes and three points collinear are not considered). 
Ihis procedure is then used to create the triangular elements which, when assem- 
bled. form the hat function (see Figure 5). For a given grid point. say (4,7), six of 
these triangular planes are joined. corresponding to the six triangles surrounding the 
point. NNE. ENE, SE. SSW.WSW, NW (see Figure 4). The unknown function ts 
then interpolated over the six elements. 

Once the unknown function is represented by the combination of these six 
interpolations, we can use Maple to compute the volume integrals. and the surface 
integrals of the gradients, in Equation ITI.17. The results of these integrals provide 
Lhe coefficients which comprise the FVE stencils. 

Thus, we have the following stencils for the volume and surface integrals re- 


spectively for interior grid points: 


13 327-5 16745 
fy,,(L Tra? rad - aa; Sel) Wee ey ioe, (III.18) 
167 —5 32745 
and 

23 
f,,, 2X VTrades)As = . 2j—-1 —87 2+] | Tes, (111.19) 

2) 
where the 27 resulting from integration in the y direction has been factored out. 


(Reeall that in cylindrical coordinates, 


I d= a es rdzdrdy.) 


Note that control volumes increase with radial distance from the origin, which ts 
reflected in a radial bias in the stencils. More specifically, r; = jh is the radial 
distance from the origin, with 7 the column index for the unknown matrix, / the 
meshsize. Thus, stencil values increase with distance from the origin. 

These stencils are applied to the unknown at a point on the grid, designated 


by its row/column position (k,7). Stencil entries indicate the values to be applied; 
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blank entries indicate a value of zero. Entry position indicates to which grid point 
ihe value is applied: left/right and up/down in the stencil correspond to neighboring 
potits in those directions on the grid. That is, if stencil values are replaced by the 


positions to which they apply. we have 


(A+) (A4+1,741) 
(A.J — 1) (A. 9) (A. + 1) 
Cee I tf) 
Clius. for example. the value that appears in the (A + 1, 7) stencil position is applied 
to the grid point that occupies that same position. 
Using techniques similar to those outlined in the two-dimensional example, 
boundary point stencils have been computed for the volume and surface integrals 


Meecha nic: Ur thys Grey, i ==, & == ill 


| 9o l 
he h 
384 lmaen) To.0; and — —3 2? Too + | a(r)rar; (IIT.20) 
Ai = 0 
i) i 
a ie qe Ti.0: (III.21) 
384 Zou k0, an 8 —-6 4 k,05 
6 l 
ZH Glee Rete —ie Oe 
327-5 16745 4) 
Ibe 


h 
247-8 W2j+5 87+3 | Tj, and 7 | 23-1 —87 23+! To,3- 


(D2 
For points adjacent to the far boundaries, the stencils in Equations IJI.18, III.19, I1.21, 
and [f{.22 are applied. However, since the homogeneous Dirichlet conditions dictate 
that these boundary values are zero, the resulting contribution after stencil appli- 
cation remains zero. Thus, in effect, far-boundary values do not contribute to the 


stencils for points adjacent to these boundaries. 
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We can now combine all of the above stencils to generate the operator L". 
l 7 , 
1a cen a ar a (ite?) 


where / is the step size on the grid. The imatrix representation of L*. M. is a block 


tridiagonal (NV + 1)?x(N + 1)? matrix of the form 


fe ae wee 0 0 
ie Ai oe 0 0 
Oo ae, ae es 0 0 
NiO 70 oy 0 ; 
0 
(ak en 0 i a aa 
1 o Ponies htt ah tee 0 L A | 


where A, A, and [ are (N +1)x(N +1) generic banded matrices (not all tdentical); A 
is tridiagonal, A is upper bidiagonal, and [ is lower bidiagonal. When M multiplies 
the matrix of (N + 1)? values of T, arranged lexicographically as a column vector, 
the matrices A, A, and [° produce the effect of the stencils “reaching” function val- 
ues respectively on the current row, the row above, and the row below. Numerical 
experiments using /atlab to construct the matrix M for N = 8,12,16,20, and 24, 
with g = and a = } and J, indicate that it has full rank. Thus, we expect a unique 


solution to the linear algebra problem which is a discretization of the heat equation. 
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IV. RELAXATION SCHEMES 


The FW method. with weighted-average time-stepping, las been used to 


discretize the continuum problem. 


OT 
— =KV°T. iy ) 
dt 
Ly 
a gars) i (ake (IV.2) 
Mie itl matrix form, 
M(T"*"| = f(T"), (1V.3) 
on a rid of meshsize h = y. The input for these equations is the solution at the 


current time step, T” (where T” and T” are used interchangeably), the result of 
solving the equations is the solution at the next time step, T”*+!. As was indicated 
in Chapter ITI, there is good reason to believe that the matrix M is of full rank and, 
thus. expect a unique solution to exist for the linear algebra problem that arises from 
the discretization of the continuum problem. Solution by direct methods requires 
factorization of M and, since M is both large and sparse, solution by direct methods 
mav be impractical. We therefore consider iterative methods to solve the matrix 
equation at each time step. These methods generate, for each time step, a sequence 
of approximate solutions, ent the choice of relaxation scheme determines whether 
or not the sequence {T5)} converges. Upon implementation of a relaxation scheme 
that produces a converging sequence of approximate solutions at each time step, we 
begin with an initial temperature distribution, 7a and the solutions are stepped in 
time. The desired result is a sequence of solutions, corresponding to a sequence of 


time steps, which tends to steady state. This process allows for the evaluation of 


'Here the subscript (s) indexes the sequence of approximate solutions; elsewhere, subscripts 
without parentheses indicate grid position or vector components. 


rail 


various Values for the time step. as well as of various weightings used in averaging the 


current and subsequent time steps. 


A. [ITERATION MATRICES 


It is common in constructing iterative methods to propose a splitting for the 
inatrix M = E—F. where linear systems of the form Ez = b are “easy” to solve (see 
Ref. 9}). The sequence of approximations, (iiaee is generated by ET oe = IF, ah 
aud. therefore, it is natural to define an iteration matrix. P = E7'F. so that 
‘ior = PT, E-'f. Additionally, if T* is the exact solution so that MT" = if 
then ET* = FT* + f and T* = E-'FT* + E-'f. Hence the solution. T*. is a fixed 
point of the iteration Tie = EFT.) eet 

The matrix P = E7'F is also called the error propagation matrix, since if 
ie, = PT;.) - Eas or ee = PT,,) te. (B a constant vector), and €(s41) 1s the 


error at the (s + 1)st step, then 


den = ey = ie 


—_— 


[PT;.) “5 B) — (PT* + B| (since / = Paes “ B) 


PT,., — PT 


eine te) 
and therefore 

€(s+1) = ee: (IV .4) 
Using induction, it is easy to show that €(,) = P*é). By Equation IV.4, €4) = Pé() 


and 


Ea) = Pet 
P[Pé()] 


= P*E(o). 


Zo 


— 


iv iuduction, assume that €,) = P*éo), in order to show that in) = Pere oy. 


Now. 


C(k+1) a Pea 


Bee ean) (by the induction liypothesis } 


hk = 
= Pp! mas 


Thus. since €a41) = Plt" é(), we conclude €(,) = P*é(), for s any integer. .\ddi- 
tionally. 
Esl] = WE*Zoyl] < PHM Zol (IV.5) 
which will be useful later. 
One of the sunplest relaxation methods is the Jacobi method. also referred 
(o as simultaneous displacement (for a more detailed account of all of these methods. 
see [Ref. 3] or [Ref. 9]). It ts produced by solving the /th equation of Equation 1V.3 


for the /th unknown, 7), / = 1: (N + 1)*. Before proceeding to a discussion of the 


iteration matrix, we present the component form for this iteration scheme: 


anh old 
eres) = a es eS aa (Ev" ) a eh )] 


| [V.6) 
kJ Sub ? ( 
aay a ier o)) 


Where )>y and >°y are. respectively, the sum of volume stencil entries and surface 
stencil entries applied to the current values of the unknown, daa except at the grid 


pomt on which the stencil ts centered: 


a Coie Tle eeO a as 
aoe old , (old rs 
Dy =| +027 - 178 +(327 41 T, |, (IV.7) 
+(16j — 5)Ty2y4-1 +(325 + 5) Thy 


and 
(old) 
~ old k+1,3 
O ee old 
des =] 4-14 27)T GS TWEE OME a (IV.8) 
+2577). 
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One iteration sweep consists of computing Equation [V.6 for each component of the 
unknown vector 7’. Provided that the process converges, the sweeps continue until a 
desired level of convergence ts reached. 

Another. more succinct. way to present this method is in matrix form. If 
we write the operator matrix as the sum of a diagonal matrix (D), a strictly upper 


iriangular matrix (U), and a stictly lower triangular matrix (L), 
M=D+U+L, 


ile matrix equation to be solved becomes 


—_ 


(D+U+L)(T] =f. 


Now define Ey = D and Fy = —(U + L), so that this may be written as 


D[T| BYUiEE IEF <2 jf en 


E,(T] = F,(T| + f 


OT as 


31 
| 


=De (OL) i ee eemor 


aa) 
| 


E;'Fy(T] +Ej;'f, 


which corresponds to solving the (th equation for JT), / = 1:(N + 1)*. If we define 


the Jacobi iteration matrix as 


P, = E7'Fy, Or 
= Ds (Uren 


the matrrx form of the Jacobi method becomes 


feet) = pee ae (IV.9) 
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\ modified version of this method. called the weighted Jacobi method, ts 
determined by introducing a weight. U< 2 ~ 1 (w= 1 is the original Jacobi method ) 


(see (Ref. 3]). The matrix form is 
P(new) = Pi es we a (IV.10) 
where [is the indentity matrix. E, = D. and F, = (| —w)D —.u(U +L), and 


a> cli 


= D“'{(l-—w)D —u(U + L)} 


(| —w)I+wP, 


is the weighted Jacobi iteration matrix. In this method, the new approximation Is a 
weighted average of an intermediate approximation and the old approximation: the 
ltermediate approximation is the Jacobi iterate of the old approximation. 

The weighted Jacobi method computes all of the components of the new ap- 
proximation before it begins to use them in the next approximation. This requires 
storing both the current and new approximations; a simple modification allows for 
updating the current approximation “in place’, using only the storage required for 
one approximation. The Gauss-Seidel method incorporates the new changes as soon 
as they are computed by overwriting the current approximation with the new (see 
(Ref. 3]). More important. however, is that (on model problems) the Gauss-Seidel 
method generally converges about twice as fast at the Jacobi method (see [Ref. 9]). 
In component form, this iteration scheme Is 


~ Laisi(Ev) - 2h (0s) 
9945 ee anh) 





Tae Ih aes (IV.11) 


are 


where the arrow indicates overwriting, and 


‘ old : old) 
(325 — 5) TE, +165 + 5) Ta 
|) O27 yy FeO Ot th eats eV at2) 
+(16j —5)Th") | 4-827 +) TS? 
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and 
iden 
Dee =| H-1 +27) H1 +2 |- (Vas 
ONE 
In matrix form. using M = D+U0+45L. Eg =(D+1L), and Fq = —U. we have 


— 


(D+L)[T] = -U[T]+f, or 
Fo(T] +f, 


les 
au 
wel 

[| 


Te) (DEL TUTC} +4 (D+4L)'f, or 
Tinew) te Ba Hoes oe Bae 


This corresponds to solving the (th equation for 7), using the new approximations for 


components 1,2,....¢— 1. Now, define the Gauss-Seidel iteration matrix, 
Pe: Ea Fe 
= —(D Ar Eiew 


so that the iteration scheme in matrix form becomes 
je) Pe eV | me yo ie (IV.14) 


With the above lexicographic ordering of the (N + 1)? components of T and the 
components updated in ascending order, the effect of a Gauss-Seidel sweep ts to start 
at r = 0 and update in the radial direction for each vertical step, starting at z = 0 
(see Figure 8). 

As with the weighted Jacobi method, we can make a modification to the Gauss- 
Seidel method. Define a parameter 7 € R (7 = 1 is the original Gauss-Seidel method) 
and the modified method is successive over-relaxation, SOR, (see [Ref. 9]) which, 


in matrix form, is given by 


Trew) — pT!) 4 vB! fF, (IV.15) 
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ay ja 5 O Cnn. 
kK 8 8 FF One) Fe 
kl sg sg £4 a 6hOCUtltéKs 





: © old approximation 


Mi onew approximation 


Figure 8. Gauss-Seidel sweep. 
where, with E, = D+ yL and F, = (1 —7)D—yJU, we have 


P, 


elo) De!) 


De 


Similar to the weighted Jacobi method, SOR is a weighted average of an intermediate 
approximation and the old approximation. 

The question of interest now ts whether or not the sequence of approximations. 
ey}. generated by Bee, = FT.) + f, converges. Therefore, we make use of the 


following theorem: 


Theorem [V.1 Suppose that if ER" andM =E-F € R**" 1s fg If 
E is nonsingular and the spectral radius of E~'F satisfies p(E~'F) < 1, then the 
iterates Te , defined by ET, (st1) = = FT.) ae converge to T = M- Lf for any starting 
vector Vey 


The proof is found in [Ref. 9], and makes use of the following lemma: 


Lemma IV.1 /f p(E7'F) < 1, then (E7'F)* + 0 as k 4 00. 
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Table I. Spectral Radii for Crank-Nicholson Time Stepping (a = $). 


| N=3|N =120N = 16 i S20 
| Jacobi | .9499 | .9709 | .9801 
= 9817 | = 9925 | 9943 | .9955 | 
=i et | oe ee 





Jt 


Table I]. Spectral Radi for Implicit Time Stepping (a = 1). 


The matrices M, Ej, and Eg, and iteration matrices, P, for various values 
of N, a = 3 and 1, and g = h, have been constructed using Matlab. We have 
experimentally verified that M, Ej, and Eg are nonsingular, and the spectral radii of 
the error propagation matrices, P. have been computed. Due to memory limitations 
with Watlab, the calculations are made for relatively small values of N. The results 
are indicated in Tables I and II, where the value of @ indicates the type of time 
stepping, and the values of w and y are the weightings in the weighted Jacobi and 
SOR methods respectively. 

The results of these numerical experiments indicate that both p(P,) < 1 and 
p(Pa) <1, with p(Pg) < p(Pz). The modifications to the Jacobi and Gauss-Seidel 
methods do not provide any improvement; the spectral radii for both methods are 


higher for the modified iteration matrices than for the original iteration matrices. 
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Lliis. the Gauss-Seidel method. used with either the Crank-Nicholson or implicit 
Hime-stepping scheme. appears to be the best chotee of these relaxation schemes. [How- 
ever. if the spectral radins of P is near unity, convergence inay by unacceptably slow 
since the error tends to 0 like p(P)*. as indicated by the lemma and [lel] < !! PIS!) Coy] 
(equation 1V.5). As is evident from the above tables. even for a moderately spaced 
ne eV = 10. Goan ad lanS , or 1), p(P) > .9, and the spectral radius 


icreases with .V. 


B. ALTERNATIVE SCHEMES: LINE RELAXATION 


Both the Jacobi and Gauss-Seidel methods give rise to iteration matrices. are 
implemented in a straightforward manner. and are attractive because of their sim- 
plicity. While the ease of implementation is a significant advantage, the convergence 
properties of either or both may not be acceptable and, therefore, alternative schemes 
may be desirable. One such, which seems reasonable based on the geometry of the 
problem, is line relaxation. There are two obvious options in this regard: radial or 
vertical line relaxation (see Figures 9 and 10), where either an entire row or column 
is updated simultaneously. Both options require the solution of an (V + 1)x(.V + 1) 
tridiagonal matrix for each row/column update in the unknown matrix. That is. while 
the previously outlined relaxation schemes (point relaxation) proceed by successively 
solving a collection of algebraic equations for one unknown, line relaxation proceeds 
bv solving a succession of matrix equations. For example, suppose that the matrices 
A, A. and [ are tridiagonal, upper bidiagonal. and lower bidiagonal, respectively, and 
that T; and f; are the portions of the respective unknown and right hand side vectors 


in Equation IV.3 that are rows in their corresponding matrices: 
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Nacdial line relaxation consists of solving the following succession of matrix equations: 
TY = An = 25). 
T, = A (f,—1T, Ai) 
Ui 1 SNE), ain! 
Ty = AT!( fy — TS). 


[nu order to illustrate how line relaxation works relative to the FVE stencils, 


let 
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so that. for the radial and vertical relaxations respectively, we must solve 
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simultaneously for each row or column in the unknown inatrix. 
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Figure 9, Updating an entire row at one time. 


The computational cost of updating an entire row/column at a time is higher 
than a row/column update using either the Jacobi or Gauss-Seidel methods. However, 
this type of relaxation may be sufficiently efficacious to warrant the extra expense, in 
that convergence may be achieved significantly faster. For comparison, we now con- 


sider local mode analysis of the Jacobi and Gauss-Seidel methods and line relaxation. 
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Figure 10. Updating an entire column at a time. 


C. LOCAL MODE ANALYSIS 


While an examination of the spectral radii of iteration (error progagation) 
matrices can be instructive, it is often useful to conduct a more detailed analysis. 
What follows 1s a local mode analysis of various iteration schemes. Since the error 
propagation matrices indicate how the error evolves during the iteration process, we 
use a DF'T (see Appendix A, [Ref. 10]) to expand components of the error equations 
associated with the various relaxation schemes. The coefficients in these expansions 
are the factors by which the corresponding modes of the error are magnified /reduced 
for each relaxation sweep. Thus, by examining these coefficients, we can determine 
how quickly the error tends to zero, which indicates how quickly the solution con- 
verges. The following analysis does not apply to points on the boundary, it is only 
valid for interior points. Thus, while it gives more information about the functioning 


of the scheme in the interior of the domain, boundary peculiarities are not addressed. 
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We begin by recalling that the current error. €(,), is the difference between 


— 


the exact solution, 7". and the current approximation. 1(5), 


where we desire that the sequence. {€,)}. tend to zero. Substituting this expression 
into Equations [V.6 and [V.11, we have the following error equations for the Jacobi 
and Causs-Setdel methods where. in order to avotd confusion with the exponential 


fuuction in the DFT expansion. we represent the components of € by o,.,: 
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where we a Se and ae are now applied to the vz,, (see Equations IV.7.1V.8.1V.12, 
and [V.13). Equations [V.22 and IV.23 indicate how the sequence {€,)} evolves 
during the iteration process. [f we expand these relationships in a discrete Fourter 
transform, (DFT) (see Appendix A or [Ref. 10]), the magnitude of the transform co- 
efficients will indicate the “amount” of each mode of the error that is present in each 
component of €(,). We then compare coefficients to determine the ratio between the 
new and current approximations for each component of the error. In other words, the 
DET allows us to analyze the growth of the error by examining component behavior. 


Expanding Equation I[V.22 in a DFT, where 
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and the superscripts (0) and (n) indicate the components of the old and new approx- 


lmations, we have 
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Making use of the orthogonality property of the complex exponential (see Appendix A) 


we can equate individual terms of the sum, and then divide by C(I, m) to give 
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[1 order to determine the greatest factor by which a mode of the error 1s multtplied, 


we seek the maximum of 





nae over the values /,m = -+ +1: *. In other words, we 
seek to determine a mak on how well we can expect this type of relaxation scheme 
to perform. This ratio is a function of N,d,m and 7 and, while difficult to determine 
analytically, may be calculated numerically. Using A/atlab, the maximum of this ratio 
lor (1 = 0: N has heen calculated for N = 16,32; « = 3,1; and 7 = |, * N—1; the 
results are indicated in Table III (see page A for a discussion of the equivalence 


of centered indices, 1,m = —% Be ae les and non-centered indices, /,m = 0: .V). 


~ ? 
Additionally, by considering the matrix of grid points /,m = 0: N, we can determine 
a correspondence between sample points on the grid and oe, - rh frequency 


(see Figure !1 and Appendix A). The matrices of values a 





v a | for j= N-—1 are 
depicted in Figures 12 and 13. 


The information in Table III indicates that the Jacobi method results in a 
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-(n) 
L | . . 
maximum of ao ereater than one: this value does not appear to depend on grid 


Larr 





position. In other words. there is at least some component of the error that is magni- 
lied. instead of reduced. by this iteration scherne. Moreover, it appears that the value 
inereases with the number of intervals on the grid. Additionally. Figures 12 and 13 
iudicate that the maximum of this ratio occurs for both low and high frequencies. 


Thus. it appears that the Jacobi method will not be effective in generating a solution 


to the point-source problem. 
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Figure 11. A two-dimensional representation of frequencies corresponding to a single 
set of non-centered sample points. 





In order to compute the ratio we for the Gauss-Seidel relaxation scheme, we 
rewrite Equation [V.23 as - 
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Fieure 12. Ratio of amplitudes of new to old Fourier coefficients, Jacobi relaxation 


10 4 
5 
a=), Y= . j = N —1, Frequency ranges are as in Figure 11. 


mifere 4 = 16)—5, B = 327-45, C = 32)—11.D = 224), F = 32j94-11.F = 32j;-5,6C = 
pees, and Ho= 2),) =—1-2),J = —$7,h = 14+ 29,0 = 27. Expanding ina 


DFT. equating terms in the (double) sum by virtue of the orthogonality property of 
the complex exponential (see Appendix A), and dividing by C(d,m), we have 
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igure 13. Ratio of amplitudes of new to old Fourier coefficients, Jacobi relaxation: 
NV = 32,a=1, 737 = N —1. Frequency ranges are as in Figure 11. 
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As above. we now seek to maximize | vie . over the values !,m = -= ae + in order 


to determine a bound on how well nts scheme performs. We again use Matlab to 
compute the maximum of this ratio for l,m = 0: N for N = 16,32; a= = 5, l; and 
we ~. N — 1; the results are indicated in Table IV. Additionally, the matrices of 
values of i! for 7 = N —1 are depicted in Figures 14 and 15. 

The informer een in Table IV indicates that the Gauss-Seidel method results 
in a maximum of rr less than one. Thus, all frequency components of the error 


are reduced by this Feianee scheme, which is what we seek. However, the value 
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Table [V. Maximum Values of ue for d.in = 0: N for the Gauss-Seidel Method of 


th el 


IvelaN ation. 


increases with grid size N, and also as the time stepping is shifted from Crank- 
Nicholson to fully implicit. Additionally. while all of the values are less than one, 
as we move to grids with larger .V and toward a fully implicit time scheme. they 
lhecome close to unity. This means that. although we may reasonably expect to 
see this relaxation process converge, 1t may be quite slow. Another point is that 
the maximum value depends on grid position; the shorter the radial distance, the 
larger the maximum. which apparently is a reflection of the radial bias of the stencils. 
Additionally, Figures 14 and 15 indicate that the maximum of this ratio occurs only 
over the low frequencies. and that the values associated with error components over 
the high frequencies are quite low. This performance over the high frequencies will 
be of importance in Chapter VI. Thus, it appears that the Gauss-Seidel method will 
be effective in generating an iterative solution to the point-source problem, albeit 
potentially slow to converge. 

We now apply similar techniques to the radial /vertical line relaxation schemes. 


The error equations come from rewriting Equations IV.20 and IV.21 respectively as 
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igure 14. Ratio of amplitudes of new to old Fourier coefficients, Gauss-Seidel relax- 
ation: N = 32, a= 4, 7 = N—1. Frequency ranges are as in Figure 11. 
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As above, we expand these relations ina DFT, make use of the orthogonality property 


to equate terms in the sums, and divide by C(l,m) to obtain 
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lieure 15. Ratio of amplitudes of new to old Fourier coefficients, Gauss-Seidel relax- 
oon. .¥ = 32,a¢= 1, ) = iN —1. Frequency ranges are as in Figure !1. 
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Table V. Maximum Values of a for d,rn = 0: N for Radial Line Relaxation. 
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Vhius. the numerators of the ratio meen are given by 
lm 
he :27l 12r(l-+m) aKkh :2mt 
— [— fi Ba = Cen we — —[-Le™® | 
g384 Z 
atid 
h? 327m 12m(l+m) akh 127m 
5384 Ke N —Ge Sa ee om 
and the denominators are given by | 
h3 :27(l+m} 
(ayaa Te emaa = ee eich 
aKkh —:2rl 127m , i2nm 
5 He wy +fe~r +/+ Neal 
and 
h? r _ r2mr(i+m) = a mae 127 
‘Taga ite Sai eran” are ~ +D+Fen | 
aKkh —127l 127m 120 
as [He“N +le vn +J/+Le% ]. 
vm) | 
As above, we now seek to maximize vel over the values /,m = —= N+ be x. We 


continue to use Matlab to compute the maximum of these ratios for t a = () oy ier 
N= 16532; a. 02> panda) ame wWN — 1; the results are indicated in Tables V 
and Vf. Additionally, the matrices of values of vee for 7 = N —1 are depicted in 
F rea res tiOema ys hs acanlenlioe 

The information in naples V and VI indicates that both line relaxation meth- 


lV, 
ods result in a maximum of ry et less than one; radial line relaxation seems to promise 
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Fable VI. Maximum Values of vi for ,mz = 0; V for Vertical Line Relaxation. 
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better performance. Many of the results for the line relaxation schemes parallel those 
noted for the Gauss-Seidel method. All frequency components of the error are reduced 
by the line relaxation schemes; their performance would appear to be slightly better 
than that for the Gauss-Seidel method. Additionally, the maximum values increase 
with both the grid size NV. and as the time stepping ts shifted from Crank-Nicholson 
to fully implicit. As with the Gauss-Seidel method, while all of the values are less 
than one, as we move to grids with larger NV and toward a fully implicit time scheme, 
they become close to unity. This means that, although we may reasonably expect 
to see this relaxation process converge, it may be quite slow, and that the increased 
amount of computational work as compared with work for the Gauss-Seidel method 
may not be worth the payoff. The maximum values for these schemes also depend on 
grid position. For the horizontal line relaxation, the shorter the radial distance, the 
smaller the maximum; for the vertical line relaxation, the shorter the radial distance, 
the larger the maximum. Once again, apparently, this is a reflection of the radial bias 
of the stencils. Additionally, Figures 16, 17, 18, and 19 indicate that the maximum 
of these ratios occurs only over low frequencies, and that the values associated with 
error components over the high frequencies are quite low. In fact, the values for the 
line relaxation schemes over the high frequencies appear to be about one-half the cor- 
responding values for the Gauss-Seidel method. As noted earlier, this performance 
over the high frequencies will be of importance in Chapter VI. 


Some general trends are evident from the results of the numerical experiments. 
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igure 16. Ratio of amplitudes of new to old Fourier coefficients, radial line relaxation: 
oe - ) = N —1. Frequency ranges are as in Figure 11. 


Since the maxima should all be less than one to guarantee convergence, it appears that 
the Jacobi method will not be useful in iteratively solving the point-source problem 
(Equation II.2). Of the schemes that have maxima less than one, those that have the 
smallest maxima are those that reduce error components most quickly and, therefore, 
are the schemes that will converge most quickly. The line relaxation schemes have the 
smallest maxima; the maximum values over the high frequencies are about one-half 
the corresponding values for the Gauss-Seidel method, but there is little increase in 
performance over the low frequencies. To determine whether or the extra work in 
using line relaxation is worth the effort, a cost comparison of getting to the same 
level of error as other relaxation schemes will need to be done. Additionally, the 
inaximum values vary depending on the time stepping scheme and the grid position; 
they are higher for implicit time-stepping than for Crank-Nicholson, indicating that 
convergence will take longer for the former as compared to the latter, and the maxi- 


mum values increase with N. Hence, as we move to a grid with more nodes, we not 
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ligure 17. Ratio of amplitudes of new to old Fourier coefficients. radial line relaxation: 
V=32,a=1,7 = N —1. Frequency ranges are as in Figure 11. 


only have more work to do by virtue of the greater number of equations to solve, we 
also have to work harder to reduce each component of the error. That is, the effect 
of moving to a finer grid is to more than quadruple the work required to solve the 
problem. 

The Gauss-Seidel and line relaxation schemes appear to be the most efficient 
in terms of reducing error components. Additionally, in Chapter VI. the performance 
of these schemes over the high frequencies will be of interest; we will want to use those 
schemes that have the smallest maximuin value over the high frequencies. In particu- 
lar, for the Gauss-Seidel and line relaxation schemes, the high frequency components 
of the error are eliminated much more quickly than the low frequency components. 
Considering the amount of work to be done and the performance of the scheme, we 


implement the Gauss-Seidel method in the solution process. 
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igure 18. Ratio of amplitudes of new to old Fourier coefficients, vertical line relax- 
OMe = oon 7 ) = N—1. Frequency ranges are as in Figure 11. 





Figure 19. Ratio of amplitudes of new to old Fourier coefficients, vertical line relax- 
ation: N = 32,a=1, 7 = N—1. Frequency ranges are as in Figure 11. 
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V. ITERATIVE SOLUTION 


Ihe results of Chapter [V indicate that the Gauss-Seidel method is the sim- 
plest. of the relaxation schemes that appear to result in a converging sequence of 
approximate solutions, and therefore Gauss-Seidel is our method of choice. Having 
chosen a relaxation scheme, we must now, as part of the tterative solution process. 
make suitable choices for the type of time-stepping, a, and the time step size. g; the 
goal is to maximize accuracy and minimize computational work. We now conduct 
experiments to solve Equation [V.3 by iterating with a solver which is based on the 
iteration matrix, Pg = -(D + L)7'U, and use this process to evaluate some of the 
possible choices for a and g. The technique is to solve first at a single time step by 
Herating until the difference between successive approximations, eee - ie or, 
perhaps better, ||M ae ~ f(T))|\, is negligible. This approximate solution, T”, is 
used as input for the time-stepping scheme so that 


7 N-+1)? 


- 
(6 ie 2 - f atav + - aja | Voy -ndo|T, 


(as in Chapter IIT) becomes the new right hand side in Equation IV.3. The iteration 
process 1s repeated to obtain the approximation at the new time step. The solution 
is stepped in time in this fashion until the difference between solutions at successive 
time steps is negligible, representing a steady solution. 


An algorithm for this process might look something lke the following: 
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Given an initial guess. T°. and tolerances db, d», 
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2) Sete: si ame 
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s+l 
S ese 
Eud while 
+) Tt! — ie 
. 
Return 7"*! as the solution at t= (n+ 1)At 
=n 1 
End while 
Tisteady 23 “ie 


Stop 


Starting with an initial temperature distribution of T = 0 everywhere, bound- 
ary conditions specified by Equations IJ.18, 11.19, and II.20, N = 16, and a = S 
several values of g are used in an attempt to achieve a reasonable solution. The re- 
sults of these experiments indicate that for larger values of g, say g = fh, the solution 
exhibits a non-physical oscillation in time (see Figure 20) at the origin, where the 
heat source is located. Instead. the solution should monotonically increase until it 
levels olf at steady state. 

The oscillation can be removed by making the temporal step size smaller, say 
gy = h°® (see Figure 21). This step size, however, results in negative temperatures 
along the z = 0 axis near the origin. Negative temperatures, like the oscillations, are 
physically meaningless. A valid solution process must eliminate both of these non- 
physical characteristics from the solution. However, it turns out that for a = 5 there 
is no value of g that produces a physically meaningful solution. In particular, for 
6 hz, there are both oscillations and negative values in the approximate solution 
(see Figure 22 for a comparison of several values of g). Thus, in order to compute a 


realistic solution, it is necessary to use values of a > >: 
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Figure 20. The maximum values of the solution for each time step, where a = 5 and 


(pa 

In an attempt to produce a solution that is physically meaningful. numerical 
experiments are conducted using values of <a <1, and values of g ranging from 
h? to h for each value of a. Solutions with no oscillation and no negative values are 
possible for, say a = 3 but this requires that the temporal step size be as small as 
—— hsiz. The size of the temporal step indicates how much work will have to be 
done to reach a steady solution: the smaller the step, the longer to reach steady state. 
[In order to minimize the work, we must maximize the temporal step size. Further 
experimentation indicates that for a = | (fully implicit time stepping) and g = A, 
the result is not only a physically meaningful solution (1.e., no oscillations and no 
negative values), but also a reasonable temporal step size. Therefore, the iterative 
solution is computed using these values. 

Iterative solutions for N = 8,16, and 32, are computed using the same initial 
temperature distribution and boundary conditions as above. One measure of accu- 


racy of these solutions is determined by comparing them with an analytic solution 
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Figure 21. The maximum values of the solution for each time step, where a = . and 


=a 

(Equations I[.15 and II.16). However, the problem that gives rise to the analytic so- 
lution presented in Chapter II] 1s somewhat different from the linear algebra problem 
that we are solving. Therefore, as in |Ref. 11] we approximate the exact solution of 
[;quation [1.2 by solving Equation IV.3 for N = 64 (see Figure 23). We then compare 
this solution with solutions of Equation IV.3 on coarser grids (N = 32 and N = 16) 
to get approximations of the discretization errors, which may be used to give an in- 
dication of the order of the accuracy of the solution. The measure that we use is the 


discrete energy norm, defined by 
|D* || = (LADS, D*)s, (V.1) 


where (-,-) denotes the Euclidean inner product, L" is the operator defined in Equa- 
tion [V.2, and D* = T" — T* approximates the discretization error, where T” is the 
solution to Equation IV.2 on grid h, and T* is the “exact” solution (i.e., solution of 


Equation IV.2 with N = 64 sampled on grid h). (see [Ref. 11]) 
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Figure 22. The niaximum values of the solution for each time step, where a = O 
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indicates g = h®. 
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indicates g = h, “+” indicates g = h?. indicates g = h?, 
Care must be taken in analyzing the figure since the apparent equivalence of time 
steps cau be misleading, e.g., for g = h* and h = <, it takes 16 time steps to “equal” 
one time step for g = h. 


The solutions for an mnital time step', for a time of one second, and for steady 
state have been compared using the discrete energy norm, with the resulting dis- 
cretization errors indicated in Table VII. We might hope to see a common factor of 


decrease as we move to finer grids. That is, if the discretization error on grid /f were 


h 


2 


of the order e = ch?, for some constant c, then the ratio of errors for grids h and 


'With g = h, the size of a time step differs with grid size. We have adjusted for tls difference by 
requiring that the same “amount” of time elapse for solutions used in the comparison, e.g., for the 
initial time step, the “exact” solution on N = 64 after eight time steps is compared to the solution 
on N = 32 after four time steps, to the solution on N = 16 after two time steps, and to the solution 
on N =8 after one time step. 
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Figure 23. The steady state “exact” solution for N = 64. 
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Thus, a decrease of a factor of four would indicate that the process achieves O(h?) 
discretization error (see [Ref. 11]). Table VII indicates that the discretization error 
decreases by a factor of about two when inoving from grid size N = 8 to N = 16, 
and by a factor of about four when moving from grid N = 16 to N = 32. A grid 
size N = 8 may be too small to adequately reflect the accuracy of the solution, which 
would leave the factor of four, possibly suggesting that the discretization error is 
O(h*). Even so, a specific estimate of the discretization error is best deferred until 


more information can be obtained, i.e., solutions on finer grids. 
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Table VII. Discretization Errors for the Iterative Solution. 





Vale MULTILEVEL SOLUTION 


Now that we can solve a discrete representation of the heat equation (lqua- 
tion IV.2) iteratively, we want to apply multigrid techniques in order to attempt to 
accelerate the solution process. Multigrid is a method to improve on solution by 
relaxation by making use of the advantages of working on successively coarser grids 
(for a more complete treatment. see [Ref. 3]). It has been used with marked sue- 
cess In solving a variety of problems, specifically elliptic partial differential equations. 
Although the problem we are solving is parabolic, in the time-stepping regime we 
inust solve an elliptic problem at each time step. Thus, multigrid may prove useful 
in streamlining our solution process. We begin by introducing the multigrid method 
and some of the exigencies of the use of multiple grids. We then analyze some of the 
characteristics of the method using local mode analysis, and present numerical re- 
sults from the solution process. The amount of work to achieve the iterative solution 
serves as a baseline against which the computational work for the multilevel solution 
Is compared. 

As indicated in the analysis of the Gauss-Seidel and line relaxation schemes 1n 
Chapter IV, the high frequency (oscillatory) components of the error are extirpated 
much more efficiently by these schemes than are the low frequency (smooth) ceompo- 
nents. The result of applying a relaxation scheme to generate an approximate solution 
on a specific grid size (call it a fine grid, size N), is that the oscillatory components 
of the error are smoothed. After sufficient work has been done on the fine grid to 
smooth the error (in effect, when the relaxation process stalls), the problem may be 
shifted to a coarser grid (grid size *) where the smooth components of the error 
beconie oscillatory (see figure 24, taken from [Ref. 3]). The relaxation scheme is then 
applhed again to smooth the oscillatory components of the error. The information 


gained from smoothing on the coarse grid is transferred back to the fine grid, where it 


becomes a correction to the original approximation. That is, the result of smoothing 
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ou the fine grid, transferring to the coarse grid. smoothing on the coarse grid, and 
transferring back to the fine grid 1s to produce a coarse grid correction. This pro- 
cess of using multiple grids to obtain an approximate solution has the advantages of 
a) more effectively targeting the error components. and b) requiring less work, since 
the coarse grid has only 2~¢ the number of points on the fine grid, where d is the 
dimension of the domain. Multigrid also allows for ways to improve upon the initial 
vuess that is used in the smoothing (see [Ref. 3]). However, since we are solving a 
time-stepping problem, where the current solution becomes the initial guess for the 
hext tune step, this will not be a concern for us. Moreover, while there is a good deal 
more to multigrid than coarse grid correction, this concept forms the foundation of 


our solution process. 


k = 4 wave on N = 16 


k =4 wave onN =8 


Figure 24. A wave with wave number k = 4 on 2" (N = 16) is projected onto 07 
(N = 8S). Porgy = 1674) 4 implies tinainthe ae 1S 7 = the way up the spectrum, 


while, for N = 8, k = 4 is the wave that is 5 = 7 way up the spectrum. Thus, the 


coarse grid “sees” a wave which is more oscillatory. 


In order to make effective use of this technique, we must specify what informa- 


tion to transfer, and how to transfer it; first we consider what information to transfer. 
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Reeall that we are attempting to solve Iquation [V.3. 


and suppose that we have an approximate solution ©; the error is given by € = [7 — 6. 


— 
ar 


where £* is the exact solution. Thus. the errer satisftes 


Me = f—Mo 


III 


fe Cy) 


\ 


—_~ 


where ris the residual. The smooth components of the error are the troublesome 
ones, sinee the relaxation schemes “kill” the oscillatory modes. Since smooth modes 
ona fine grid become oscillatory when projected onto a coarse grid, it 1s natural to 
consider transferring a representation of the error, 1.e.. the residual, from one grid to 
the next’. In this way, the relaxation schemes used on the coarse grid will reduce 
components of the error that could not be reduced on the fine grid. Additionally, we 
know that relaxation smooths the error, and therefore we can accuratcly represent 
the error on the coarse grid. After transferring the residual to the coarse grid. we 
can relax ou the coarse grid version of the residual equation (Equation VTI.1), solving 
for the error. Thus, when we have deterntined an approximation to the error on the 
coarse grid, we can transfer it back to the fine grid as a correction to the current 
approximation. That is, stnce T° =0+6, if we know @ we can correct it by adding 
an approximation of €. This process can be outlined as in the following steps, where 
h represents the grid spacing on the fine grid, H the spacing on the coarse grid, and 
we assume 2h = A. 


Coarse grid correction (two-grid scheme) 


e Relax on M"T? = f" on 2" to obtain an approximation v". 


'There are a number of other good reasons to transfer the residual, as outlined in [Ref. 3]. 
However, there are other multigrid schemes that do not transfer the residual, such as the Full 
Approximation Scheme (FAS) (see [Ref. 2] and [Ref. 1]). FAS is useful in dealing with nonlinear 
problerns but, since our problem is linear, we will not employ it. 


63 


e Compute the residual 7 = f? — M*v?, 


e Solve the residual equation M4%e7 = F? on Q” to obtain an approximation 


to the error e!. 


e Correct the approximation obtained on * with the error estimate obtained 

on DN: o* — G 4+ e” (Ref. 3I 
The task of solving on 2” can, thus, be exchanged for the task of relaxing on 
Q” and then solving on 2%. The requirement to solve on 2/7 can be treated in like 
fashion: the task of solving can be exchanged for the task of relaxing, and then solving 
on the next coarser grid. This procdeure of transferring to successively coarser grids 
can be continued until a grid 1s reached on which an “exact” solution is possible. 
That is. the solution can be generated by recursive use of the coarse grid correction, 


which can be represented as follows: 


Solve on 2” by relaxing on 2’ and solving on 7", 
Solve on 27" by relaxing on 2?" and solving on 94", 
Solve on 2%" by relaxing on 04" and solving on 92". 


“Exact” solve on coarsest grid. 


Correct on 4". 
Correct on 22". 
Correct on 12". 


This process of starting on a fine grid, moving to the coarsest grid, and then going 
back to the fine grid is known as a V-cycle (see Figure 25). Each time the transfer 
1o a coarser grid is made, the relaxation process there smooths another portion of 
the error component spectrum (see Figure 26). Thus, by the time the coarsest grid 
is reached and the problem has been solved, all of the error components have been 
smoothed. When the problem is solved on the coarsest grid, the approximate error 
is transferred to the next finer grid to become a correction. The corrected error 1s 


then passed to the next finer grid, and so on until we are back on the finest grid, 
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where the original approximation 1s updated with the composite correction from the 
coarser grids. While there are other niultigrid methods which may prove more useful. 


depending on the nature of the problem (see [Ref. 3]), we use the V-evele. 
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Figure 26. Error component spectrum for various grid sizes. On the fine grid. hf, 
relaxation smooths the high frequency components of the error (heavy black line). 
When the problem ts transferred to the next coarser grid, 2A. the high frequencies 
of the remaining (unsmoothed) components (heavy black line) are smoothed. This 
process continues all the way down to the coarsest grid, where the problem ts solved, 
eliminating the rematning components of the error. By continuing to transfer to 
successively coarser grids, all frequency components of the error are smoothed. 


So far. we have not indicated how information will transferred. nor how an 


error estimate on the coarse grid (solving M"%é" = Fr) will be computed. We will 
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deal with the former question in the next section; the latter question revolves around 
deciding what to use as the coarse grid operator, M”. One approach is to simply 
discretize the problem on the coarse grid in precisely the same fashion that it is 
discretized on the fine grid. This has the advantages that the work has already been 
done, and that it 1s easy to implement computationally. While this method generally 
works (on model problems), it has the disadvantage that it 1s often not true to the 
physics, e.g., conservation may no longer be enforced. There are other methods for 
determining the coarse grid operator which are true to the physics and allow for strong 
theoretical treatment of convergence and other properties (see [Ref. 2] and [Ref. 3]). 
However, because discretization on the coarse grid generally works. and because it 
simplifies the computations. forming M” in this manner is the method that we use. 

We conclude this section with a discussion of what happens to boundary con- 
ditions when we transfer the residual. The boundary conditions for our problem are 
a mixture of Dirichlet and Neumann conditions (Equations II.18, II.19, and I1.20). 
By relaxing on Equation IV.2, 

UN = 


(where we have dropped the superscript /) we obtain an approximation v, which gives 
rise to the residual equation, r = f — Lv]. Suppose that T* is the exact solution, so 
that f = L[T*]; we require that the approximation satisfy the same conditions as the 


exact solution on the boundary, 02Q. Thus, the residual equation becomes 


r = LT] — Lv 


] er ae 
= JG ~ anv?)T-av) LG V?)vdV] 


= - (T* —v)dV — ak fowr — V*v)dV 
Gg JV 
| Bs 3 oO: 
= =a —v)dV —an f (VT Vv) ads. (V1.2) 





Therefore, since the boundary conditions satisfy (T* = v)|aq and (2 = 2)|aq’, all 


22 = V-n denotes the outward normal on @2. 
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boundary conditions for r become homogeneous; the first term on the riglit: hand stde 
of Equation VI.2 indicates that this ts true for the Dirichlet conditions, the second 
term indicates that it 1s trne for the Neumann conditions. Thus, since we are solving 
the residual equation on all the successively coarser grids. we impose homogeneous 


boundary conditions on all but the finest grid. 


A. INTERGRID TRANSFERS 


ln order to transfer inforination between grids. we develop operators that 


restrict data from a fine grid to a coarse grid, and interpolate data from a coarse 


SUPerscripl 
subscript 


ord to a fine grid (see [Ref. 3]). These operators are designated by I , where 
the subseript indieates the grid from which the information is being transferred. and 
the superseript madieates the ertd to which the information ts being transferred. Thus. 
the restriction operator ts indicated by If’, since the information goes from the line 
erid to the coarse grid, and the interpolation operator ts represented by Ij), since 
the information goes from the coarse grid to the fine grid. In other words, restriction 
is a process of determining what values to assign to coarse grid points, based on the 
values at fine grid points; interpolation is a process of determining what values to 
assion to fine grid points, based on the values at coarse grid points. In [Ref. 1], such 
operators are developed which take advantage of FVE characteristics. We eonsider 
Iwo types of restriction operators, one which follows from the physics and one which 
is very simple, and one type of interpolation operator. 

One of the advantages of the F VEE method is its fidelity to conservation consid- 
erations. The restriction technique that we describe first follows from the conservation 
notion. In order to determine what value to assign to a coarse grid point, the fine grid 
is laid over the coarse grid, where the coarse grid control volumes align with fine grid 
lines (see Figure 27). The control volume for a given coarse grid point includes all or 


part of the control volumes of nine fine grid points. Hence, the value of a quantity on 


the coarse grid includes contributions from the value of the quantity on each of the 
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nine fine grid points. For each of these fine grid points, the contribution 1s determined 
by computing the fraction of the associated fine grid volume that is contained in the 
coarse grid volume. The value of the quantity at the fine grid point is multiplied by 
this fraction, and the result is the contribution to the coarse-grid value. The value 
of the quantity at the coarse grid point 1s thus the sum of the nine contributions 
determined in this manner. 


Fine Grid 





grid lines 


Be oan eee control volumes 


Coarse Grid 





Figure 27. The “conservation restriction” operator. 


In the case of cylindrical coordinates, care must be taken when deciding what 
percentage of the fine grid volumes fall within each coarse grid volume. Since the 
volumes change with the radius, the restriction operator will be a function of grid 
location. The stencil for the restriction operator gives the percentage of the fine grid 
value that is assigned to the coarse grid value based on the percentage of the fine 
grid volume that falls within the coarse grid volume (see Figure 28). For example, 


the amount of the value of the point at (4,7 +1) that is used in the sum is found by 
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Figure 28. The fine grid region, comprising contributions from the nine fine grid 
control volumes, that corresponds to a coarse grid control volume. The central fine 


orid point is labeled (A, 7), (k,j; even) for which the radial distance r = jh, and 
corresponds to coarse grid point i =k 


lirst determining how much of its fine grid control volume (right-center, or rc) falls 
within the region that corresponds to the coarse grid control volume. This volume 
is then divided by the total control volume for the pomt (4,7 + 1) to determine the 
volume percentage. So, tf the fine grid location is at the radial distance r = jh, then 


we have for the volume re (a toroidal prism) 


Or 


Or 
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Similarly. the ratio for the left-center (lc) portion of the stencil is 


and the remainder of the stencil is given by 
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This is a stencil for interior points; similar calculations determine the stencils for 


boundary points. 
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Figure 29. Neighboring coarse grid points. 


Since this operator is both somewhat out of the ordinary and dependent on grid 


location, we consider it worthwhile to try to get a sense of how this restriction operator 
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treats the quantities being restricted. For tlus purpose. we consider a comparison 
between the values assigned to neighboring coarse grid pomnts, as in Figure 29, when 
Lhe restriction Operator acts on a grid with constant value 7;, = 1. The value for 


coarse grid pomt mmmber | is given by 
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Thus, if the restriction operator acts on a constant vector, the values assigned to the 
coarse erid pomts decrease as the radial distance increases. 

Does this make sense in terms of the physics of the problem? Should a constant 
function remain a constant after restriction, or not? In general, one might expect 


that restriction would transfer a constant to a constant. The conservaton restrictlou 


el 


operator acts on the basis of volume percentages and, thus, considers the values 
assigned to each control volume in terms of the actual volume. In some sense, this 
restriction operator converts the values assigned to a control volume into a “density”, 
and then transfers this amount. In the case of cylindrical coordinates, assigning a 
constant value to each control volume means that the ratio of the value to its volume 
decreases with radial distance. Thus, the result of the above computation is no 
surprise. Moreover, it seems to make sense physically to transfer “densities” in terms 
of enforeing conservation. However, in order to regain the fine-grid interpretation of 
the information after it has been transferred may require more work. In particular, 
some compensation may need to be made relative to the coarse grid volumes (..e., 
convert the “density” back to the original type of information). The question of how 
or Whether to make this compensation will not be addressed here. Additionally, it is 
certainly possible that the information to be transferred does not make sense in the 
context of densities. In that case, physical intuition may need to be suspended while 
the efficacy of the solution process is determined. 

The other restriction operator we introduce 1s the injection operator. Injection 
is accomplished by simply assigning values to the coarse grid points directly from the 
corresponding fine grid points, ug; = v},9; (see [Ref. 3]). The injection operator 
has the advantage of being a linear operator, whereas the conservation restriction 
operator may not be linear. Note that, when the conservation restriction operator 1s 
applied to a grid with constant value 7;,,; = 1, the result is 72,9; = 4. If the injection 
restriction operator 1s applied to the same grid, 74.2; = 1. Thus, care must be taken 
when comparing these two restriction operators. 

In order to transfer from the coarse grid to the fine grid, we must choose an 
interpolation operator. Since we are assuming that the basis functions for the FVE 
method are linear, a natural choice is linear interpolation. A continuous H-piecewise 
linear function is also a continuous h-piecewise linear function. That is, the “coarse” 


finite element space is a subspace of the “fine” finite element space (see [Ref. 1]). 


te 


This means that linear miterpolation correspouds to the “uatural embedding” from 
the coarse fintte element space to the fine element space. Its advantages are its 


sunplicity and its linearity. The stencil ts given by 
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B. AMPLIFICATION MATRIX 


Now that we have our fine grid operators. coarse grid operators, and intergrid 
transfer operators, we conduct analysis using the DFT (see Appendix A) similar to 
that in Chapter [IV to obtain some indication as to how a two-level multigrid scheme 
using these operators will perform. We start by noting that the process of relaxing 
on the fine grid, obtaining a coarse grid correetion, and relaxing again on the fine 
grid can be represented by the two-level error propagation matrix, P’?(CG)P". Note 
that, while h as a superscript identifies grid spacing, 4; and v2 are exponents. Thus, 
P’ and P'? represent relaxation v; and vy times on the fine grid. and CG 1s the 
coarse grid correction. Recall that coarse grid correction 1s the following series of 
Steps: 


Coarse grid correction (two-grid scheme) 


e Relax on L*T* = f* on * to obtain an approximation v". 


e Compute the residual r* = f? — L?v*. 


e Restrict the residual rf = /f'r*. 

e Solve the residual equation L%e =r on QO” to obtain an approximation to 
the error e”7, 

* Interpolate the error e* = [/'e". 


e Correct the approximation obtained on 2" with the error estimate obtained 
eee areas aces 


iC 
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where / is the identity operator. J? and /;, are the restriction and interpolation 


' is the coarse grid smoother. Each 


operators, L’ is the fine grid operator, and (L””)~ 
of the operators in the two-level error propagation matrix is expanded in a DFT 
(see Appendix A and Chapter IV), making the resulting expansions functions of 
Pa —< +1: ~. the indices of the expansion in the frequency domain. These 
eXpansions are used to construct an amplification matrix, A(/,m), which is also a 
function of / and m. The largest spectral radius of this matrix over the values that 
ime have on the coarse erideae 1.72) — = +1: (this relationship is discussed 
later). will give us a two-level asymptotic convergence factor, \. This factor, 
much lke the values determined for the relaxation schemes, will give a bound on 
low well the two-level scheme can be expected to perform. In order to guarantee 
convergence, \ < | is necessary; the smaller \ is, the better. 

In Chapter IV, we used DF'T expansions of the error equations derived from the 
relaxation operators to determine the frequency domain coefficients for the operators. 
We now use the same technique, again using error equations, to determine the matrix 


entries for the operators in the two-level scheme. The amplification matrix 1s formed 


by multiplying together these operator matrices, 
A(!,m) = (P*)?(1- 1, (17) 1 LP"), (VI.3) 


where 


id 


e L” is the 4x4 matrix for the fine grid discretization operator. 

e P* is the 4x4 matrix for the fine grid relaxation operator. 

* (L?)-} is the 1x1 matrix for the coarse grid relaxation operator. 
. iy is the 1x4 matrix for the restriction operator. 


’ ie is the 4x1 matrix for the interpolation operator. 
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e Tis the 4x4 identity matrix. (see [Ref. 2]) 


We go through the detains of computing the matrix entries for the fine grid operator, 
L”: results for the remaining operators are then presented. 


From the FVE stencils for L”. we have the operator equation 
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+4) V94-1.2, 
where the superscript fA indicates a vector component on the fine grid. (old) and 
(new) indicate before and after the operator acts, and subscripts 2A and 2) are the 
respective fine grid row and column indices. These indicial values are used because 
we shortly will need to discuss the correspondence between the fine and coarse grid 
points; only those fine grid points with indices 24h and 2) correspond to coarse grid 
points (with indices 7 and /:). Ordinarily (see [Ref. 2]), the consideration of indicial 
values is not a concern. In this case, however, since the stencils are grid location 
dependent, and since the analysis of the amplification matrix must apply to both 
the fine and coarse grids, the 7 that is used is the column index on the coarse grid, 


requiring that 27 be the column index on the fine grid. Thus, expanding both sides 
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where 2 < 27 < NV —2 on a grid labeled 0: N, since we are in the 1teriong 
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When we transfer information from the fine grid to the coarse grid, we must 
remember that not all of the freqencies on the fine grid can be represented on the 
coarse grid. The highest frequency that can be resolved on any grid is the Nyquist 
frequency (see [Ref. 12]), f = 54, where Az is the grid spacing. Hence, on the 
eoatse @rids Ween * and f = *. Therefore, in order to make this analysis 
germane to both fine and coarse grids, we must rewrite the sums in the expansions 


so that |l], |m| < * This can be done as follows: 
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Phus. rewriting the DFT expansion of the operator equation to consider those 
lrequenees that can be represented on the coarse grid aud, making use of the or- 
theeouality property of the complex exponential (see Appendix A), we eam equate 


individual terms of the sun, and then divide by C(/, im) 4 0 to give 
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The matrix of coefficients for the fine grid operator is a 4x4 diagonal matrix, 
L' = (A,;) (see [Ref. 2]). To see why the matrix is 4x4, recall that the sum in the 


original DFT expansion of the operator equation was rewritten so that the indices 





satisfy |l|,|m| < *. One result of this is to rewrite each of the individual terms of 
the sum as a combination four terms (Equation VI.4). That is, for each dimension 
of the problem, represented by / and m, there are two components of each individual 
term in the sum, which are indexed by |/|, |m| << 2 and 2 < ||, |m| < ¥. Thus, the 
matrix of coefficients is a 2?x2* matrix, where each row in the matrix, say the zth row, 
comprises the coefficients used in determining the corresponding zth component in the 
new vector as a linear combination of components of the old vector. In this case the 


inatrix 1s diagonal, since for each of the four components of the new vector, the only 


contribution from the old vector comes from the corresponding component. In other 


18 


Jy\re She 
words. the only contribution to f ‘ , comes from Nac 


1) and VO are the 4x] vectors ieee components have indices (f.7). (f+ 4 


(dom + *), and (/ + ~. m ++). The four diagonal entries 





he inl 
sea (61) — 53 4. (32) 
eal Fe (O87 — je Be Yo) asm ha 


+-(64j — 1l)e7e™ 





ani 











AS) eee Pete = 
akKkh t27l 127m 
——— (4je™ + (-1 +4 je 


\ 13 ( (64) 5) 227i (32 E 4 5) 
Aggy = |—{—-(647—-—5)e% -(- 5)e 





+(647 — Il)e 
12m(i+m) 


+448) + (647 + Ll)e* 


| 


=16j (1 + Aye + 4je )] 


So, (() = LA 
a mo ; 


mee eer ee 


er 


1:27 rn 





1:2 7(l+m} 


“x + 448; + (64) + 11)e oe 


ea ee ale 





akh e or 


s— (dye + (-1 + 4y)er™ 


ol 


227i 





—167 + (1 +4j)eon" —4djeN ) ; 





\ a ((64) Je —(32j +5)e 
= =e 
cae g354 / J 


—(64j — jes 





eae 











—(327 — 5)e7 + (64) + 5)e7 
/ F :27m™ 
——— (je = (—1 4+ 4j)e7 


2 





—167-—(1 + Qe x at 4je—™ )| 


and 


27m 





+ 448) — (327 4 11)ex 


ent) 





” 


where 


V 
2 


ick a 


h . - +27l Z r2m(l+7n) 
en a sh —5)e™ +(327+5)e 
g384 
—(64j —ll)e"* 4+ 4485 — (647 + l1)eoe 


+(32) —5)e7 ae — (647 + 5)e 4) 


aKh = 


= (—4je* —(—1 page 


216) = eee “xe — dye “w Y). 











A similar process must be applied to determine each of the other matrices 


that comprise A(l,m). In order to compute the entries for the matrix P”, recall from 


Chapter IV, Equation IV.25, that we can write the error equation, e”*” = Pe", as 
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similar calculations to those above, we have the four entries, I];;, of the diagonal 


tnatrix P?: 
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In similar fashion, the entry for the Ix] matrix (L!)7!, corresponding to the 


coarse grid relaxation operator, is given by 
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The equation for the conservation restriction operator, [}!, is 
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tou) = aoa) Cet Ca) and expand this relation in a DFT. Considering 


those frequencies that can be represented on the coarse grid and, making use of the 


orthogonality property, we can equate individual terms of the sum. We divide by 


= w2rki = we2mpmn . 3 
(‘(l,m) =e wn e-N_, the coarse grid analog to C(l,m), to give the compouents of 





the matrix of coefficients corresponding to the restriction operator, [7 = (A,;). This 
matrix is 1x4 since the output of the operator, a single component on the coarse grid, 
is a linear combination of 2? components on the fine grid. The entries in the matrix 


are as follows: 
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The equation for the injection restriction operator is Ui = Came which, using 
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The linear interpolation operator gives rise to the following relations: 
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[expanding these in a DFT, considering the coarse grid frequencies, and making use 


of the orthogonality property, we can equate terms of the sums and divide by C(I, m) 
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The components, 0;), of the matrix of coefficients, Le. are found by solving the above 
system of four equations for the fine grid components, such that Vi = 1 ee The 
matrix I, is 4x1 since the output of the operator, the 2? components of the fine grid 


vector, is the result of interpolation operator acting on the single component of the 
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coarse grid vector. That 1s, 
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Thus, we have all the elements to construct the amplification matrix, A(1, 7). 


(equation V1.3). We now have the tools to analyze the expected perfomance of a 


two-level scheme. 


C. NUMERICAL RESULTS 


We now conduct a number of numerical experiments to investigate the perfor- 
mance of the multigrid method, again using Matlab. First we examine the behavior 
of the amplification matrix, A(/,7m), and the largest spectral radius of the family of 
matrices A(l,m) for all l,m. We then present the results of the repeated use of the 
V-cycle to solve Equation IV.2, as compared with the iterative solutions. 

[n order to make use of the amphfication matrix, we must choose a relaxation 
method, and the number of times to relax on the way down and on the way back up. 
i.e., We must choose 14, and v2 in Equation VI.3. We experiment first with Gauss- 
Seidel relaxation and three different combinations of (1,12) V-cycles: a (0.1) cycle, 
a (2,1) cycle, and a (10,10) cycle (recall that a = 1, meaning fully implicit time 
stepping, and the time step g = /). The first two types of cycles are fairly common 
(see [Ref. 3]), while the latter is included for reference. The information presented 11 
Tables VIII, IX, and X indicates the asymptotic convergence factor, A, for specific grid 
sizes, grid positions, and type of restriction operator, either conservation or injection. 
Linear interpolation is used as the interpolation operator in all of these experiments. 
Tables XI and XII indicate the asymptotic convergence factors for horizontal and 
vertical line relaxation and a (2,1) cycle. 


Table VIII indicates the asymptotic convergence factors that result from a 
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(0,1) cycle, Gauss-Seidel 
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Fable VITL Asymptotic Convergence Factors for the (0,1) Cycle Using Gauss-Seidel 
Relaxation. 


single iteration sweep on the fine grid coupled with the course grid correction for 
both the conservation and injection restriction operators. In general, the conservation 
restriction operator seems to indicate better performance, since its convergence factors 
are lower than those for the injection operator (recall that the asymptotic convergence 
factor is the largest spectral radius of the family of amplification matrices, which 
comes from the DFT expansion of each component in the two-level error propagation 
matrix). Similar to the analysis on relaxation schemes in Chapter IV, the convergence 
factors decrease with radial distance, indicating an improvement in performance, but 
Increase with grid size. Even though the convergence factors for all three grid sizes are 
less than one, indicating that convergence is very likely (but not necessarily certain, 
since these factors tell us nothing about what is going on at the boundaries), at least 
some of the factors for all three grids are relatively high, indicating that convergence 
may be quite slow with the (0,1) cycle. A point of interest is to compare these 
results with those in Table IV, which indicate the maximum ratio of new to old error 
components for Gauss-Seidel relaxation. One might expect the results in Table VIII 
to be universally better than those in Table IV, however this is not the case. In 
fact. on a grid of size N = 32, the largest asymptotic convergence factors for both 
the conservation and the injection restriction operators are larger than the largest 
ratio for just a single Gauss-Seidel sweep on the fine grid. In other words, a single 


Gauss-Seidel relaxation sweep is predicted to do better than a (0,1) cycle. This is a 
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Table IX. Asymptotic Convergence Factors for the (2,1) Cycle Using Gauss-Seidel 
Relaxation. 


fact, worth remembering when the multilevel solution is analyzed. 


Table LX indicates the asymptotic convergence factors for a (2,1) cycle; as 
expected, the factors are lower than the corresonding factors for the (0,1) cycle, and 
are universally lower than the ratios associated with a single Gauss-Seidel sweep as 
indicated in Table IV. Ilowever, another comparison is to ratse the ratios in Table 1V 
to the third power, corresponding to three relaxation sweeps on the fine grid. In other 
words, since for a (2,1) cycle there are two relaxation sweeps on the fine grid before 
the course grid correction, and one afterward, for a total of three relaxation sweeps 
on the fine grid, it seems reasonable to compare the performance of a (2,1) cycle to 
three iterations of the Gauss-Seidel method, represented by raising the Gauss-Seidel 
factor to the third power. For example, in Table IV for N = 32. 7 = 1, the ratio 
is .9567 which, when raised to the third power ts .8756. This value is lower than 
either of the values in Table IX for N = 32, 27 = 2. This also holds true for the 
asymptotic convergence factor associated with N = 32 and 2) = : (this plenomenon 
does not occur for N = 32, 27 = N —2, nor for N = 16). Thus, the overall predicted 
performance of three Gauss-Seidel relaxation sweeps for N = 32 is better than the 
predicted performance of a Gauss-Seidel (2,1) cycle. 

Similar to the information in Table VIII, the performance predictions in ‘Ta- 


ble IX improve with radial distance, and worsen with increase in grid size. Also, as 
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Table X. Asymptotic Convergence Factors for the (10,10) Cycle Using Gauss-Seidel 
Relaxation. 

above, the conservation restriction appears to predict better performance than injec- 
tion. The factors are still relatively high, especially for the NV = 64 grid size, so that 
convergence will likely be slow. 

The (10,10) cycle is included for reference so that the performance of the 
two-level cycle can analyzed using a large number of iteration sweeps on the fine grid. 
The characteristics of the (10,10) cycle are similiar to those of the (0,1) cycle and the 
(2.1) cycle; the performance improves with radial distance and worsens with increase 
in grid size. The asymptotic convergence factors are smaller than for previous cycles, 
but this is expected since a good deal more work is represented by the (10, 10) cycle. 
However, using a comparison similar to that in the discussion of Table 1X, in Table [V 
for N = 32, ) = 1, the ratio is .9567 which, when raised to the 20th power is .4126. 
Here. as in the case of a (2,1) cycle, this value is lower than either of the values 

1 Table X for N = 32, ») = 1. This holds true for all the factors associated with 

= 32, but is not true for N = 16. Thus, as before, the predicted performance 

of Gauss-Seidel relaxation for N = 32 is better than the predicted performance of a 
Crauss-Seidel two-level cycle. 

Tables XI and XII indicate the asymptotic convergence factors for horizontal 
and vertical line relaxation (2,1) cycles. Similar to the results of the analysis on the 
horizontal and line relaxation schemes in Chapter IV, both line relaxation schemes for 


the (2,1) cycle indicate better performance than that for the Gauss-Seidel method; 


| (2.1) cycle, Horizontal Line Relaxation 
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Table NT. Asymptotic Convergence Factors for the (2.1) Cycle Using Horizontal Line 
Relaxation. 
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Table XII. Asymptotic Convergence Factors for the (2,1) Cycle Using Vertical Line 
Relaxation. 


horizontal line relaxation is the better of the two. One difference is that in Tables V’ 
and VI, the predictions for horizontal line relaxation worsen, aud the predictions for 
vertical line relaxation improve, with radial distance, whereas in Tables XI and NII 
the predietions for both line relaxation schemes improve with radial distance. Using 
calculations similar to those in the discussions of Tables IX and XX, the values in 
Tables V and VI are raised to the third power and compared with the factors in 
Tables XI and XII. The result is that the predicted performance of three sweeps of 
the horizontal line relaxation is better than that for the (2,1) cycle for both V = 16 
and N = 32. The predicted performance of three sweeps of the vertical line relaxation 
is better than that for the (2,1) cycle only for N = 32. These results are consistent 
with those of the Gauss-Seidel relaxation for all of the (1,12) cycles tested. That 
is, for N = 32, the predicted performance of ordinary iteration is better than the 


predictions for any of the corresponding two-level scliemes considered here. 
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In the initial stages of using V-cycles to solve Equation IV.2 at a single time 
step, we experiment with the use of conservation restriction and linear interpolation. 
The result is that this combination does not result in a coarse grid correction. That 
is. the effect of applying the two-level scheme is to generate an approximation that 
is not as good as that obtained by relaxation on the fine grid alone. Therefore, we 
experiment with the combination of the injection restriction operator and the linear 
interpolation operator, with the result that this combination does produce a coarse 
erid correction. A correction scheme using V-cycles with these intergrid transfer 
operators is implemented to solve at a single time step, and then the solution 1s 
stepped in time until a “steady” solution was obtained. 

Multigrid V-cycle solutions to Equation [V.2 are computed on grids N = 8, 16, 
and 32, using (2,1) cycles, Gauss-Seidel relaxation, a = 1 (fully implicit time step- 
ping), and time step g = h, where the initial temperature distribution and boundary 
conditions are the same as in the iterative solution. We again use the discrete energy 


norm (Equation V.1) to measure the accuracy of these solutions, 
|LD°|I| = (LDS, DY), 


where (-,-) denotes the Euclidean inner product, L” is the operator defined in Equa- 
tion 1V.2, and D*® = T” — T* approximates the discretization error, where T” is the 
solution to Equation [V.2 on grid h, and T* is the “exact” solution (i.e., solution of 
[Equation IV.2 with N = 64 sampled on grid h) (see [Ref. 11]). The results presented 
in Table XIII are almost identical to those obtained for the iterative solution (see 
Chapter V). 

As another measure of how well the iterative solutions match the multilevel 


solutions, the discrete L? norm of the solution differences, d”, 
jh = I d* d” 5 
I< | = A y ) ’ 


is computed for each of the above grid sizes and time value. The results are presented 


in Table XIV; while the norm of the differences is generally small, it increases with 
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Table XTTE. Diseretization Errors for the Multilevel Sebution. 





Rable XIV. Norins of Differenees Between Iterative and Multilevel Solutions. 


erid size. This might suggest a normalization problem, however, the factor of 7 in the 
diserete L? norm is used specifically to avoid this type of problem. More investigation 
is required to determine why there is not better agreement between the iterative and 
multilevel solutions. 

Figures 30 and 31 indicate the number of iterations per time step for the it- 
crative solutions and the number of V-cycles per time step for the multigrid solution 
for N = 32. [t is tmportant to note that, while both processes require a decreasing 
amount of work for each time step out to about 150 steps, the multilevel process re- 
quires about 5 more time steps than the iterative process to reach a steady solution: 
the multilevel process requires one V-cycle per time step for time step > 150. In 
addition, Figures 32 and 33 illustrate how the norm of the difference between suc- 
cessive solutions behaves as the solution is stepped in time. As would be hoped, the 
difference between successive solutions decreases with time stepping. However, as 1s 
evident from these latter figures, the time stepping process stalls somewhat before 
the steady solution is reached (the tolerance used is machine ¢€ = 2.22/ — 16). While 
the number of time steps to reach a steady solution on grids NV = 8,16 (not shown) 


is about the same for both the iterative and multilevel solutions, on grid Wo = 32, 


CA 
cap) 


as already indicated, the number of time steps needed for the multilevel process to 
reach steady solution is about : longer than that required for the iterative solution. 


I-ven so, the two processes on all three grid sizes stall at about the same time step. 


Iterations per time step 
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Figure 30. The number of iterations per time step for N = 32; 170 time steps needed 
to reach a steady solution. 


We conclude our analysis of the multilevel solution with a discussion of how 
much computational work must be done to reach a solution. For this purpose, we 
pick a time step, for both methods, for which the time stepping process has stalled. 
For NV = 8, use time step s = 50; for N = 16, use time step s = 80; and fone yaar ee 
use a time step of s = 150. The norm of the difference between successive solutions 
for these time steps is on the order of 107'°. We use calculations similar to those in 
(Ref. 3] in computing the cost. If we consider a work unit, WU, to be the cost of 
one relaxation sweep on the finest grid, then we can estimate how many work units 
are associated with each V-cycle. Since we are using a (2,1) cycle, relaxation occurs 
three times on each level; one relaxation sweep on the grid 2?" requires p~? of a work 


unit, where d is the dimension. Adding these costs, and using the geometric series as 
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V-cycles per time step 


Mieure 31. The number of V-cyeles per time step for V = 32; 230 time steps needed 
tu reach a steady solution. 


an upper bound, we have 
| | 3 
V-cycle Cost = 341 + Se Se ane \WU < =o. 


Thus, for the two-dimensional case we have Cost = 4WU. A comparison of the 


computational cost for the various grid sizes 1s given in Table XV. 
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Total Cost (WU) 
4984 16304 93124 
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Table XV. Computational Cost for the Iterative and Multigrid Solutions (Not Includ- 
ing the Cost of the Intergrid Transfers). 
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Norm of difference between successive time step solutions 
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igure 32. Norm of difference between successive iterative solutions for N = 32; 170 
time steps needed to reach a steady solution. 


The cost of implementing the intergrid transfers is not included in the com- 
putations presented tn Table XV. In considering the amount of work represented 
by one work unit, counting “adds” and “multiplies” each as one floating point op- 
eration (flop), the number of flops for one relaxation sweep on a grid of size N is 
approximately 41N%. Injection restriction requires no arithmetic; linear interpola- 
tlon requires approximately aN" flops. Thus, the cost of the intergrid transfers for 
the two-level scheme is about (32° /41N?)WU = =.WU. Again using the geometric 
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series as an upper bound, we have that the cost of the intergrid transfers is given by 


3 


3 
aie 42-64 .-.-42-™ WU < —WU. 
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Transfer Cost = 





Thus, we update the information in Table XV and present the results in Table XVI. 

While not dramatic, the use of multigrid does result in a savings of compu- 
tational cost. The cost savings increase with grid size; the savings on a grid of size 
N =8 is only 20%, on a grid of size N = 32 the savings is almost 50%. In consid- 


ering these results, several questions need to be answered. For example, are these 
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Norm of difference between successive time step solutions 
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Pieure 33. Norm of difference between successive multilevel solutions for NV = 32: 230 
Lime steps needed to reach a steady solution. 


Total Cost (WU) 
V-cycle cost (WU) 4992 16354 93286 


Iteration cost (WU) 25631 | 103185 





Table XVI. Computational Cost for the Iterative and Multigrid Solutions, Including 
the Cost of the Integrid Transfers. 


results consistent with the local mode analysts and the mformation that comes from 
the amplification matrices? Does the multigrid solution show “typical” performance: 
does the process work as one would hope? 

In general, the results are consistent with the local mode analysis, which 
predicts that convergence will be slow. The prediction of the amplification matri- 
ces is difficult to apply, since the asymptotic convergence factors apply to two-level 
schemes and not to V-cycles. Nevertheless, the convergence factors did not predict 


that the two-level scheme would provide improvement over the iterative process on 


93 


erid size .V = 32, and therefore these results are somewhat of a pleasant surprise. 
In other words. multigrid performed better than might have been expected based on 
the two-level convergence factors. However, typical multigrid performance reaches 
convergence (to the level of truncation error, see [Ref. 3] and [Ref. 2]) in just a 
few V-eveles. Thus, the multilevel solution developed here is disappointing, since 
convergence requires a very large number of V-cycles. 

There are no obvious causes for the disappointing performance of the multi- 
level solution of Equation H].2. One of the problems may be the use of cylindrical 
coordmates, which results tn a radial bias in the FVE stencils. The discretization of 
the time derivative using finite differences, while using the FVE method discretize 
the spatial derivatives. may also have some tmpact on the solution process. It 1s not 
clear, however, that either of these affects the multilevel solution any differently than 
it affects the iterative solution. The choice of relaxation scheme is a major factor 
in determining the convergence of the process; perhaps the use of line relaxation, or 
some other type smoother, will improve the performance. Again, however, it is not 
certain that this will result in an improvement of the multigrid performance relative 
to the iterative process. The determination of the coarse grid and intergrid transfer 
operators is also something of a concern. This is one area that clearly affects the 
performance of the multilevel solution as compared to the iterative solution. The fact 
that the predicted performance of iteration sweeps alone is better than the predicted 
performance of the two-level schemes emphasizes the need to consider the use of dif- 
ferent coarse grid and intergrid transfer operators. Finally, the point-source problem 
includes a singularity at the boundary. Most of the analysis has been done for interior 
points, and the solution results do not drastically differ from what is predicted by the 
analysis. Neverthless, the treatment of the boundaries, particularly for this problem, 


may provide a means to improve the solution process. 
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VII. CONCLUSION 


Phe purpose of this work is to apply the FVE method to discretize the axisym- 
metric heat equation, and inplement multigrid in solving it. In addition to using the 
FV method to discretize the spatial derivatives, finite differences are applied to the 
tine derivatives, resulting in a time-stepping scheme that makes use of a weighted 
average of the current and future time steps. The use of cylindrical coordinates results 
ina discretization stencil that has a radial bias, the effect of which is evident in the 
analvsis of various relaxation schemes. Gauss-Seidel is the relaxation scheme chosen 
lor implementation im the iterative and multilevel solution processes. Due to numer- 
ical anomalies encountered in computing the iterative solution, the weighted-average 
tine-stepping scheme becomes a fully implicit backward time-stepping scheme. Its 
use requires more computational work than the weighted-average scheme, but its use 
is required to achieve a physically meaningful solution. The specifics of the multilevel 
technique are then presented, including the development of the coarse grid and inter- 
erid trausfer operators, and the predictions of expected performance for the two-level 
scheme. The multilevel solution is then implemented; its results are analyzed and 
compared to the results of the iterative solution process. The results of our work are 
somewhat disappointing in that we are not yet able to achieve the hoped-for level of 
success. There are, however, a number of positive results that come from this work, 
specifically, the application of the FVE method to a problem in cylindrical coordi- 
nates and the use of Afaple software to compute the necessary integrals. Additionally, 
we know that the usual multilevel techniques do not produce the expected results and 
therefore there is ample room for further study. 

There are several possibilities for future work. One specific question that has 
not been addressed is the determination of the consistency of the discretization oper- 


ator. Numerical results are useful in attempting to verify this result. but consistency 
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should be proven analytically. In conjunction with this idea, there may be a require- 
nent to review the assumption that the solution can be adequately represented by 
a collection of continuous, precewise linear basis functions. There ts no evidence as 
vet that would require this review, but it is possible. Additionally, we are assum- 
ing a uniform step size that apphes to both the r and z directions. One possible 
modification to this approach ts to implement a non-uniform grid, perhaps one that 
results in control volumes that are the same for every point on the grid. The point of 
this modification would be to elimimate the radial bias in the discretization. For the 
discretization of the time derivatives. the requirement to produce a physically mean- 
inefiul solution has resulted in a fully implicit time-stepping scheme. While there may 
be another method to discretize the time derivatives. it will have to meet the same 
criterion and, therefore, may not be much of an improvement. 

All of the above considerations apply generally to the problem, and not specifi- 
cally to multigrid. There are at least three areas in which multigrid performance may 
be improved: intergrid transfer operators, coarse grid operator, and treatment of 
boundary conditions. The intergrid transfer operators used are the simplest available 
that result in coarse grid correction and, consequently, choosing different operators 
may improve performance. Other possibilities include the use of a restriction operator 
that is based on the piecewise linearity of the basis functions, or perhaps an inter- 
polation operator that is based on enforcing conservation. In other words, choose 
intergrid operators that satisfy the variational condition, I}, = c(I#)’, for c some 
constant (see [Ref. 3] and [Ref. 1]). The choice of the coarse grid operator may 
need to be made in conjunction with these intergrid operators. That is, despite the 
computational complexities that would be introduced, choose a coarse grid operator 
such that the Galerkin condition L” = I, LIP (see (Ref. 3] and [Ref. 1]) is sat- 
isfied. Additionally, special treatment of the boundaries may be required, since the 
boundary of the point-source problem has a singularity at the origin (see [Ref. 2]). 


Finally, the consideration of the point-source problem is the first step in ap- 
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vlying the PVE discretization and multilevel solution to the Navier-Stokes equation 
Ihat arises from the molten-pool problem. Che mtricacies of the molten-pool problem 
oresent several challenges. ticlided are the requirement to solve a svstem of three 
PDs in three unknowns and the requirement to keep track of the moving plase- 
change boundary as the inolten pool expands. Additionally, high flow velocities and 
stall local length scales result when convection is vigorous ({Ref. 13]). causing fur- 
ther numerical compheations. Tt is anticipated that the lull Approximation Scheme 
(IAS) (see (Ref. 2] and [Ref. 1]) can be used to effectively treat the non-linear nature 
of the problem. and that the Fast Adaptive Composite Grid (FAC) method (see [Ref. 


I} will be useful in overcoming difficulties that arise in the geometry of the problem. 


oF 





APPENDIX A. THE DISCRETE FOURIER 
TRANSFORM 


Below we briefly list the definition and some of the properties of the Discrete 
Fourier Transform. DFT. lor au exhaustive treatment see (Ref. 10]; what appears 


in this appendix closely follows that treatment. 
The DFT is, in at least one interpretation, a way to discretely approximate 


the Fourier transform. For exainple. assume we are working on a temporal or spatial 


: Pe : a eI. e: oe ro pe , ‘ 
domain (—4, 5] with grid spacing Aw = A/N and grid points r, = yAr, and its 


assuciated frequency domain [—*. 4 with grid spacing Aw = 27/A and erid points 


» = mw. where jm = —*+1:4. Suppose v, denotes the sampled values. v(.r,). 


of a funetion defined on the domain [—4. 4] and that d(w,,) is the Fourier transform 


of v at the frequency grid points, w,,, where the Fourier transform is given by 


(iw) = / Clale res odie 
—900 


Using the Trapezoidal Rule to approximate the integral. we have 


N 


ras Z 
- imme l mi2nym 
B(Wm ) =e v(irje 4 dx BAK — S- Ure ar ; 
-4 N v 
g=— S41 
where m = —* +1: ~. Given the set of N sample values of v;, the DFT comprises 


the \ coefficients 


9 : 


a 


N 
( —s2mjym N N 
Vn = — Sy vje NN , for m=——+4+1:— 
rf 2 
me) 


Thus, the DFT can be considered to approximate the Fourier transform o(w,) by 


(Oe = A Vi . 


We now give a formal definition of the DFT’. 


'There are several definitions of the DFT, as indicated in [Ref. 10]. 
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Definition A.1 Let N be an even positive integer and {v;} a sequence” of N complex 
numbers where j = —“%+1:%. The discrete Fourier transform of the sequence {v;} 
is another sequence of N complex numbers, {V,,}, whose elements are given by 


N 
Ver L Ss oe (A.1) 
i) V 49 , . 

; pal 
Z 


Ni 
OD ts se 


w|=S 


The DFT may be defined for N odd as well, however we will not need this definition 


for our discussion. The choice of indices, —* +l: ~ as opposed to 0: N —1, 1s 
made deliberately in order to simplify some of the analysis that is to follow from 
applying the DFT, and because it is more natural, if less familiar. While the choice 
of indices does not affect the applicability of the transform, care must be taken to 
avoid confusion over which grid points correspond to what type of frequency (see 
Figure 34). We use the operator notation D{v;} to denote the DFT of the sequence 
{v;}, and D{v;}m to indicate the mth element of the transform, D{vj}m = Vin- 

In general, the output of the DFT, {V,,}, is a complex-valued sequence. The 
interpretation of the DFT coefficients ts essentially the same as that given for the 
(continuous) Fourier transform. Also, V,, 1s even more closely related to cn, the 
Fourier series coefficient. In many ways, the DFT ts more naturally viewed as an 
approximation to ¢m than to b(w). The mth DFT coefficient V,, gives the “amount” 
of the mth mode (with a frequency w,,) that is present in the input sequence v,. In 
contrast to the use of modes of all frequencies, as in the Fourier transform, an N-point 
DET uses only N distinct modes, with roughly ~ different frequencies. The modes 


can be labeled by the frequency index m, and each mode has a value at each grid 


point 2; where 7 = _< +1: * Therefore we denote the jth component of the mth 


DEFT mode as 


as —t2mjym 
Wry =e N ; 
lee Anat — —+ +1: =. 


~The terms sequence and vector are used interchangeably to denote the input/output of the DFT. 
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kK 
N A N N 3.N 
- —— () = —* N 
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{High trequency Low trequency High frequency 


N 


Figure 34. Using centered indices (bottom figure), the low frequencies correspond to 
hi\ < ~ for a single set of sample points; the high frequencies correspond to |f:| > = 
With non-centered indices. low frequencies correspond to0 <k < * and aN ae Se oh 
for a single set; high frequencies correspond to ~. oh re aN 


The DEFT allows us to transform from the temporal (spatial) domain tnto the 
frequency domain. The Inverse Discrete Fourier Transform, IDFT, permits us 
to transform froin the frequency domain back into the temporal (spatial) domain. We 


now give a definition for the IDF T. 


Definition A.2 Let N be an even positive integer and {Vj} a sequence of N complex 
numbers where in = — + +1; ~. The inverse discrete Fourier transform of {Vin} ts 
another sequence of N complex numbers, whose clements are given by 


N 
é t27zm 


i Uk Vane @ (A.2) 


m=-241 


As with the DET. the IDFT is defined for N odd, but we will not need this defini- 
tion. Additionally, we use the operator notation D7'{V,,} to denote the IDFT of the 


sequence {V,,}, and D7'{V,,}; to indicate the jth element of the inverse transform. 


De ven }j = U5: 


LO] 


The discrete Fourier transform and its inverse have many useful properties. 
Below are listed some of the more important ones that we will be using. First, 
lhowever, we give the discrete orthogonality property for the complex exponential, 
which is essential in working with DFTs. 


Property A.1 Orthogonality. Let 1 and m be integers and let N be a positive 
integer. Then 


se a a Noén(I — im), (A.3) 
4s 


where Sv(k) is the modular Kronecker delta. defined by 


tL) — 1 ufk =0 ora multiple of N 
n(k) = 0 otherwise 


Although we will not specifically need it, we include the relationship between 


the DFT and the IDFT: 


Property A.2 Inversion. Let {v,;} be a@ sequence of N complex numbers and let 


D{v;} ={V,,} be the DFT of this sequence. Then Dally (=) ae 


The remaining properties that we will need are presented below. Proofs and 


detailed explanations are found in [Ref. 10]. 


Property A.3 Periodicity. The complex sequences {v;} and {V,,} defined by the 
N-point DFT pair (A.1) and (A.2) are N-periodic. That is, 


Vv; =vjen and Vz, = Vmen for all intergers 7 and 7m: 


Property A.4 Linearity. /f {v;} and {w;} are two compler-valued sequences of 
length N and a and 3 are two complex numbers, then 


Diav; + Bwjtm = AD{v;}m + BD{w;}m. 


Property A.5 Modulation. [f the elements of the input sequence {vj} are multi- 
plied by oe for k a fixed integer, then the DFT of the modulated sequence 1s given 
by 

D{vjwny }m = D{U;}man = Vink. 
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Property A.6 Shifting. /f the wuipul sequence {v;} es shifted (or translated) hk units 
to the right, then 

j . ae 

Die, on 7 Von Bie 

We make extensive use of the DIT in analyzing our solution process. The DIT 

is an extremely useful tool in predieting the performanee of relaxation sehemes (Chap- 
ter IV) and in evaluating various elements in the multilevel solution (Chapter VI) via 
local mode analysis. Of the properties outlined above, the orthogonality and shifting 


properties are the two that figure most prominently in local mode analysis. 
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APPENDIX B. INTERPOLATION TOOLS 


Below are the interpolation tools developed by David Canright [Ref. 8] for 


computing the FVIe stencils using \aple. The indexing in the following program 


does NOT follow the indexing in the text: the subscript ¢ is the column indicator 


and corresponds to the column indicator j used in the text: the subscript j is the 


row indicator and corresponds to the row indicator / in the text. Thus, the indices 


(2,7) inthis Maple program indicate the 7th row and the zth column, in the text the 


indices (A, 7) indicate the Ath row and jth column. 


The routine “plane” returns the function for the plane through the three given 


poiuts. There is no error checking, sv it has probleins with special cases: all 3 poimts 


collinear (Infinite solutions); and vertical planes (no solutions). 


Dl atic a= oproclc ln lez Ng 222 Konvt 2 Ota ay cox. 
unapply( 
subs( 
solve( { 
Ged lees ieee 
ax r2+b* y2+c 
iin ee oa ca Ne a eae er 
{4,5,c}), 
axxc+b*ey+c), 
a) 
end: 


| 
by 


To interpolate the indexed function 7 on the 6 triangles surrounding 


be . 


(1. *h,j * h), use the 6 functions (counter clockwise from northeast) “tri: 


the point 


ene, nne, 


uw, wsw, ssw, se”. Note: these return functions of two variables (for example r, 2). 
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triene := (1,j.2) — > plane(i7lhy Wizp. yy ee ie 
(j4+1)*h,z. ne): 

trinne := (i,j,z) — > plane(i*h,j*h,z._p,i*h,(j+1)*h,z.n,(i+1)*h, 
(j=l) h,z-ane): 

trinw :=(1j,2) — > plane(*h, jz) hj Ze ie 
J Ane 

triwsw := (i,j,2) — > plane(i*h,j*h,z._p,(i—1)*h,j*h,z.w,G—1)*h, 
(j-1)*h,z.-sw): 

trissw := (i,j.z) — > plane(i*h,j*h,z._p,i*h,(j—1)*h,z.-s,(i—1)*h, 
(j—1)*h,z.-sw): 

trise := (1,j.z) — > plane(i*h.j 7 z.20,1 lj) ee ee 
| Wegeees: 


To convert the notation from subscripted back to indexed, use “toindex(expression, 
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variable)”, where the variable is to become indexed. 


toindex := proc (expr,var) local n,m; 
subs( 
zip( (x,y) — > (x=y), 
[var. (_sw,-_s,_se,_W,-p,-e, nw,n,ne)], 
[seq (seq (seq (varli+n,j+m], n= —1..1),m= —1..1)]), 
expr ); 
end: 


Now we define shorthand for limits of integration about the general point (2,7). 
Here. rdiag means the diagonal line, where r depends on z. (This assumes integrating 
(PES tenes mune tiie: | 


rp = ih; rw = (—1/2)"hy rec= 0041/2) hoadiae jy hh Ze 
Zp i= yh; 2s *= G— 1/2) ie cee ene 


Now compute the volume integral for the unknown temperature 7’ (which is a 


surface integral of JT after factoring out the 27 from the gy integration): 


int (int (triene (1,j,T)(r,z)*r, r=rdiag..re),z=zp..zn) + 
mt (int (trinne (1,j,T)(r,z)*r, r=rp..rdiag),z=zp..zn) + 
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eCity (ie LT yt Pt Wp =z g) 4 

itt (tris (4.), ) ir. 2) aie PW, valiae 52 =75..715 |) + 

(eae watt aye Wl ee et ac tele on 27S. 71 y a 
(int ( 


Pee Tey elt e | i ieer ye ee 75 ee71))) 3 


factor(simplify(~ )); 
toindex(” ,T); 


This set of computations produces the stencil entries for the volume integral: 
: 

384 

(162 — SG ee, Ss (327 = 2D) da + (162 + 2) ae) : 





((322 + lL) e414 + (822 + 5)73 5-1 + 224075, + (322 — 11) 7,-1,, 


ln similar fashion. the surface integral of VT is computed (which is a circula- 


tion mtegral of 7 after factoring out the 27 from the » integration): 


int (+D{1] (triene(i,j,T))(re, _ re,z=zp..zn) + 
int (+D{2] (trinne(ij,T))(r.zn)*r,r=rp..re) + 
int (+D{2] (trinw(ij,T))(r,zn)*r,r=rw..rp) + 
int (-D[1] (trinw(ij,T))(rw.z)*rw,r=zp..zn) + 
int (-D[1] (triwsw(i,j,T))(rw, A rw,Z=Zs..zp) + 
int (-D[2] (trissw(,j,T))(1,zs)*r.r=rw..rp) + 
int (-D[2] (trise(i,j,T))(1,zs)*r,r=rp..re) + 

int (+D[1] (trise(,j,T))(re.z)*re,z=zs..zp); 


factor(simplify(” )); 
toindex(”,T); 


This set of computations produces the stencil entries for the surface integral: 


A 
i GW Dt det te Oe Wa a 2a ay 20 yay ) - 


fant 
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